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ABSTRACT 


V 


Fluid  turbulence  Is  of  crucial  significance  In  many  problems  of 
scientific  and  technical  Importance.  Current  developments  In  computer 
technology  offer  the  possibilities  of  solving  the  fundamental  equations 
of  turbulent  flow  In  a  vay  never  before  possible.  In  order  to  accomplish 
this  aim,  however.  It  Is  first  necessary  to  formulate  the  essential 
theoretical  concepts  In  a  suitable  manner.  This  report  summarizes  the 
progress  achieved  to  date  In  this  connection.  Various  essential  basic 
equations  are  derived,  but  the  emphasis  Is  as  much  on  fundamental  con¬ 
cepts  as  on  mathematical  details.  More  specifically,  a  method  Is 
established  for  the  computer  simulation  of  the  detailed  stationary  tur¬ 
bulence  In  a  uniform  shear  flow.  The  results  obtainable  In  this  way 
are  far  more  comprehensive  than  any  which  could  reasonably  be  obtained 
by  physical  experiment.  The  data  generated  represents  fundamental 
information  which  may  be  subsequently  analyzed  to  establish  overall 
phenomenological  characteristics  of  the  turbulence.  The  concepts  In 
this  report  should  provide  a  sound  basis  for  a  systematic,  sustained 
and  productive  research  plan.  They  have  already  been  successfully 
applied  to  a  computer  program  which  Is  now  going  Into  operation.  Results 


of  a  typical  computer  run  are  Included  and  Illustrate  qualitative  agree- 

i 

ment  with  theoretical  predictions;  It  Is  hoped  to  present  far  more  compre- 
and  definitive  numerical  results  In  future  reports. 
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SUMMARY 


A 


Pie  Problem 

For  solution  of  the  flow  In  complex  hydrodynamic  systems  at  high 
Reynolds  number  (for  example,  upward  flow  following  underwater  nuclear 
explosions),  a  basic  understanding  of  the  mechanism  of  fluid  turbulence 
is  required.  The  Reynolds  stresses  applied  to  the  mean  flow,  various 
"eddy  diffusion  coefficients  ,"as  veil  as  the  rate  of  energy  dissipation 
from  turbulence  to  heat,  all  depend  on  the  state  of  the  turbulence  Itself. 
Various  experimental  information  is  available  concerning  these  phenomena 
for  certain  particular  configurations  of  the  mean  flow,  and  various 
theories  of  turbulence  exist  which  are  applicable  to  a  few  special  cases. 
The  basic  partial  differential  equations  which  govern  all  Incompressible 
fluid  flows  (including  turbulent  ones)  have  been  known  for  many  years. 

To  date,  however,  it  has  been  impossible  to  derive  from  these  equations 
the  relevant  overall  statistical  parameters  describing  turbulence  which 
are  required  for  technical  applications. 

Findings 

In  the  current  investigation,  the  approach  to  the  problem  of  obtain¬ 
ing  this  Information  is  indirect  rather  than  direct.  That  is,  by  the 
use  of  high-speed  digital  computers  which  have  been  available  only 
recently,  the  exact,  detailed  nature  of  particular  realizations  of  tur¬ 
bulence  may  be  obtained  numerically,  by  brute-force  integration  of  the 
governing  partial  differential  equations.  This  extremely  large  amount 
of  detailed  information  concerning  the  flow  may  then  be  treated 


ii 


statistically  to  extract  the  essential  information  required.  This 
report  describes  techniques  developed  to  obtain  the  desired  detailed 
Information,  and  presents  a  few  results  from  a  computer  program  which 
embodies  these  concepts.  Methods  are  then  discussed  by  which  this  data 
(which  Is  analogous  to  a  very  large  amount  of  physical  experimental  data) 
will  be  processed. 
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1.  INTRODUCTION 


The  present  investigation  of  turbulence  originated  in  connection 
with  a  project  involving  the  numerical  simulation  of  the  mean  flow  in 
underwater  nuclear  explosions.*  It  was  realized  quite  early  in  the 
course  of  the  explosion  study  that  a  significant  fraction  of  the  energy 
released  in  the  explosion  goes  into  the  creation  of  intense  turbulence, 
and  that  the  turbulence  and  the  mean  flow  have  a  profound  mutual  influ¬ 
ence  upon  each  other.  Of  course,  it  is  possible  to  simulate  the  mean 
flow  phenomena  to  a  first  order  of  approximation,  which  is  sufficient 
for  sane  engineering  purposes,  by  ignoring  the  complexities  of  turbu¬ 
lence  and  simulating  them  rather  crudely  in  the  form  of  a  fictitious 
increase  in  apparent  viscosity.  In  fact,  some  such  procedure  is  required 
in  any  case  to  ensure  that  the  numerical  calculation  procedure  shall 
itself  be  stable. 

It  is  apparent,  however,  that  further  progress  in  the  underwater 
explosion  problem  requires  a  deeper  understand ing  of  the  associated 
turbulence.  Furthermore,  the  phenomenon  of  turbulence  has  an  enormous 
technical  and  scientific  importance  in  its  own  right,  quite  apart  from 
the  specific  application  to  underwater  explosions.  Most  fluid  flows  of 
technical  importance  are  turbulent,  so  that  progress  in  understanding 
the  fundamentals  of  turbulence  has  potentially  an  enormous  range  of 
application. 

♦Pritchett,  J.  W.,  'toACfL  -  A  Two  Dimensional  Cylindrical  Coordinate 
Incompressible  Code",  U.  S.  Naval  Radiological  Defense  Laboratory, 
USNKDL-LR-67-97,  20  Oct.  1967. 
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It  is  an  Interesting  fact  that  the  basic  partial  differential  equa¬ 
tions  which  govern  tne  detailed  motion  Involved  in  fluid  turbulence 
have  been  known  for  a  very  long  time.  However,  no  one  has  yet  succeeded 
in  deriving  from  thes'i  detailed  equations  the  resulting  overall  statis¬ 
tical  and  phenomenological  characteristics  of  the  turbulence  which  are 
important  for  technical  applications.  Typical  statistical  features  of 
these  kinds  are  mean  kinetic  energy,  mean  turbulent  stresses,  scale  of 
turbulence,  various  correlation  coefficients,  etc.  The  ultimate  mutual 
interrelationships  among  these  statistical  characteristics  of  the  tur¬ 
bulence  and  the  relevant  features  of  the  mean  flow  field  remain  shrouded 
in  mystery. 

One  reason  for  the  difficulty  is  that  the  detailed  structure  of  the 
turbulence  is  extremely  complex  in  relation  to  that  of  the  mean  flow 
and  requires  an  enormously  greater  number  of  degrass  of  freedom  for  its 
full  description.  This  fact  is  crucial  for  analysis  by  purely  numerical 
methods.  It  should  be  borne  in  mind  that  in  the  underwater  explosion 
problem,  for  example,  the  mean  flow  Itself  Is  quite  complex  and  requires 
a  very  large  number  of  degrees  of  freedom  for  its  adequate  numerical 
description.  This  is  true  even  after  advantage  has  been  taken  of  the 
polar  symmetry  which  reduce  i  the  mean-flow  problem  from  three  to  two 
dimensions.  The  mean-flow  requirements  still  tax  the  memory  capacity 
of  the  best  modern  computers.  If  ve  now  were  to  attempt  to  add  turbu¬ 
lence  effects  by  a  direct  numerical  assault,  this  would  require  a  much 
finer  mesh  and  an  extension  from  two  to  three  dimensions.  The  demands 
on  computer  memory  capacity  and  computing  time  would  be  Increased  by  a 
fantastic  amount!  Clearly,  such  an  attack  by  sheer  force  would  be 
wholly  Impractical. 

In  view  of  the  foregoing  difficulties  with  a  direct  numerical  attack, 
and  In  the  absence  of  any  demonstrably  valid  analytical  solutions,  most 
attempts  to  Incorporate  the  effects  of  turbulence  into  various  problems 
of  fluid  mechanics  have  been  based  on  a  more  or  less  empirical  approach. 
Thus  the  unknown  relation  between  the  mean  effective  turbulence  stresses 
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and  the  mean  flow  may  be  assumed  to  follow  same  more  or  less  plausible 
mathematical  form,  perhaps  one  derived  by  a  rough  analogy  with  the  known 
mechanism  which  governs  the  viscous  stresses.  Various  assumptions  of 
this  kind  can  and  have  been  postulated,  but  none  so  far  suggested  have 
been  convincingly  demonstrated  to  follow  from  the  fundamental  equations 
of  continuity  and  motion.  Consequently,  heuristic  models  of  this  type 
generally  contain  one  or  more  empirical  coefficients  whose  values  can  be 
determined  only  by  physical  experiment.  Nevertheless,  results  can  be 
obtained  In  this  way  which  do  have  a  certain  usefulness  for  particular 
purposes.  However,  such  results  tend  to  be  rather  limited  in  scope,  are 
subject  to  occasionally  serious  Inconsistencies  and  errors,  and  do  not 
provide  sufficient  insight  Into  the  real  underlying  mechanism  Involved. 

There  is,  however,  another  possible  alternative.  We  may  subdivide 
the  overall  problem  Into  several  successive  phases,  each  one  of  which 
does  lie  within  the  capabilities  of  modern  computers.  In  the  first  place, 
we  can  turn  attention  to  the  study  of  situations  In  which  the  character 
of  the  mean  flow  has  been  substantially  simplified.  All  nonessential 
complications  of  the  mean  flow  are  at  first  eliminated.  In  this  way  the 
Inherent  and  essential  mechanism  of  turbulence  is  all  the  better  Isolated 
and  displayed  for  detailed  study  and  analysis.  Secondly,  we  can  sub¬ 
divide  the  turbulence  itself  into  two  parts,  namely,  small-scale  and 
large-scale  turbulence.  The  small-scale  turbulence  can  be  studied  first, 
utilizing  the  full  capabilities  of  the  computer  to  this  end.  While  the 
numerical  data  obtained  In  this  way  Is  enormously  detailed  and  complex, 
the  significant  overall  statistical  and  phenomenological  features  of 
this  data  presumably  lend  themselves  to  summary  and  generalization  In 
sane  much  simpler  way.  Thus  the  mean  effective  turbulence  stresses 
associated  vlth  the  small-scale  turbulence  can  be  consolidated  Into  an 
appropriate  formula.  In  fact  the  analysis  given  in  a  later  section  of 
this  report  sheds  some  light  on  the  essential  nature  of  such  a  relation. 
The  resulting  formula  can  then  be  treated  as  part  of  the  known  Input  In 
the  treatment  of  the  large-scale  turbulence.  Again  the  full  power  of 
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the  computer  is  utilized  in  dealing  with  the  large-scale  turbulence. 

The  detailed  numerical  data  obtained  in  this  second  stage  can  again  be 
generalized  so  as  to  svumnarize  the  mean  effective  turbulence  stresses 
and  other  phenomenological  characteristics  of  the  overall  turbulence. 

These  results,  in  turn,  can  be  treated  as  known  inputs  in  the  treatment 
of  the  mean  flow  field,  and  can  be  applied  to  progressively  more  complex 
mean  flows. 

Of  course,  the  relationships  obtained  for  a  simplified  mean  flow  do 
not  necessarily  apply  to  mean  flows  of  more  general  types.  Nevertheless, 
if  fundamental  insight  is  our  goal,  the  natural  progression  of  study  must 
be  from  the  simpler  to  the  more  complex  cases.  Eventually,  it  should  be 
possible  to  re-introduce  some  of  tta  complicating  factors  which  are  ex¬ 
cluded  from  the  initial  studies.  The  resulting  mathematical  models 
should  gradually  become  more  pertinent  and  general.  Meanwhile,  however, 
even  the  earlier  models  can  be  used  to  advantage  as  first-order  approxi¬ 
mations  and  introduced  provisionally  into  rather  general  types  of  flow 
fields,  such  as  the  unsteady  mean  flows  involved  in  underwater  explosions, 
for  example.  This,  at  least,  is  the  philosophy  upon  which  the  present 
line  of  investigation  is  based. 

A  more  detailed  discussion  of  the  problem  of  modelling  the  mean 
effective  turbulence  stresses  is  presented  in  a  later  section  of  this 
report. 

The  principal  purpose  of  the  present  investigation  is,  therefore,  to 
begin  a  fundamental  study  of  the  basic  mechanism  of  turbulence  by  means 
of  a  numerical  solution  of  the  equations  of  continuity  and  motion  on  a 
modern  highspeed  digital  computer. 

So  far  as  the  present  report  itself  is  concerned,  Its  main  objective 
is  to  document  the  progress  that  hae  thus  far  been  achieved,  and  to 
summarize  the  key  concepts  that  have  evolved.  A  certain  amount  of  explora¬ 
tory  theoretical  work  and  preliminary  computer  calculations  were  carried 
out  in  the  early  stages  before  the  essential  ideas  were  adequately 
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developed.  It  is  not  the  purpose  of  this  report  to  dwell  on  these 
early  attempts  which  have  played  a  useful  role  but  which  are  now 
largely  superseded. 

It  is  perhaps  appropriate  to  mention  in  passing,  however,  that  some 
efforts  were  made  to  formulate  the  problem  in  wave  space,  using  concepts 
of  spectral  and  Fourier  analysis.  While  these  methods  have  certain 
attractive  theoretical  features,  it  was  finally  concluded,  nevertheless, 
that  numerical  analysis  can  be  carried  out  more  efficiently  by  a  more 
straightforward  approach  expressed  in  terms  of  ordinary  physical  space. 
However,  the  results  computed  in  this  more  direct  way  can  subsequently 
be  re-analyzed  from  the  spectral  point  of  view.  In  fact,  it  is  expected 
that  such  supplementary  calculations  of  the  spectral  type  will  definitely 
be  made  and  will  prove  useful.  (See  Appendix  D.) 

While  several  preliminary  computer  programs  have  bo  far  been  produced, 
the  first  program  which  embodies  present  guiding  concepts  to  a  sufficient 
degree  is  that  one  designated  as  "HJRBOCODE,  MARK  V."  As  of  the  date  of 
this  writing  (January  1968),  this  program  has  Just  been  developed  to  the 
point  where  serious  calculations  can  now  begin  on  a  B.oderate  scale. 

Results  of  a  typical  computer  run  are  given  later  in  this  report.  How¬ 
ever,  extended  discussion  of  actual  numerical  results  must  await  a 
future  repci t.  The  present  report  is  aimed  primarily  at  documenting  the 
theoretical  concepts  which  are  largely  embodied  in  TURBOCODE,  MARK  V. 

In  fact,  the  present  theory  goes  somewhat  beyond  that  on  which  the  above 
code  is  based.  It  is  anticipated  that  a  more  up  to  date  version  will 
ultimately  replace  MARK  V.  At  this  time,  however,  it  is  considered 
advisable  to  obtain  a  certain  amount  of  numerical  data  from  MARK  V 
before  proceeding  with  any  major  program  revision. 

The  anticipated  progress  of  the  present  research  program  can  be 
classified  into  several  successive,  but  somewhat  overlapping  stages,  as 
follows : 

(jl)  Development  of  basic  ooncepts,  equations  and  computer  program 
for  simulation  of  turbulence. 
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(2)  Generation  and  analysis  of  turbulence  data,  and  development  of 
phenomenological  models. 

(3)  Application  and  extension  of  results  to  various  mean  flow  fields 
of  technical  importance. 

The  present  report,  however,  is  largely  limited  to  the  first  item 
in  this  list.  The  second  stage  will  also  require  a  strong  analytical 
effort.  Only  then  can  the  technological  pay-off  represented  by  the 
third  stage  be  fully  achieved. 


2.  THE  BASIC  PLOW  FIE  ID 


The  simplest  possible  mean  flow  field  relating  to  the  problem  at- 
hand  is  one  which  Is  both  steady  and  1 nc empress lb le . 

The  simplest  possible  state  of  turbulence  is  one  which  is  both 
stationary  and  homogeneous.  Stationary  turbulence  is  characterized  by 
time  invariance  of  all  statistical  properties  at  all  points.  Homogeneous 
turbulence  is  characterized  by  spatial  invariance  of  all  statistical 
properties  at  all  times.  These  statements  also  imply  constant  viscosity 
over  space  and  time. 

The  turbulence  can  be  stationary  and  homogeneous  only  if  the  mean 
flow  itself  is  steady  and  homogeneous.  For  homogeneity,  the  mean  flow 
must  have  a  uniform  and  constant  vorticity  vector,  and  a  uniform  and 
constant  strain  rate  tensor .  If,  in  addition,  we  require  that  the 
streamlines  of  the  mean  flow  field  be  straight  and  parallel,  say  parallel 
to  the  x  axis,  the  mean  flow  reduces  to  the  definite  and  unique  case  of 
a  simple  uniform  shear  flow.  We  can,  without  loss  of  generality,  orient 
the  y  axis  in  the  direction  of  the  shearing  gradient.  Therefore  the 
mean  velocity  becomes  simply 

a  A 

U  -  i  [Uo  +  fy]  (2-1) 

A 

where  i  is  a  unit  vector  in  the  positive  x  direction,  and  Q  is  the  con¬ 
stant  shear  rate.  UQ  is  a  constant  which  depends  only  on  the  arbitrary 
location  of  coordinates. 

A  further  simpl if i catic-n  of  Eq.  (2-1)  would  be  to  impose  the  restric¬ 
tion  of  zero  shear  rate,  that  is 

A 

|~  ■  i  n  -  0  (2-2) 
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whereupon 


U  -  i  U  -  constant  (2-3) 

This  represents  a  uniform  mean  flow.  It  Is  in  fact  equivalent  to 
no  flow  at  all,  since  a  set  of  axes  can  now  be  defined  which  move  along 
with  the  mean  flow,  with  respect  to  vhich  axes  the  mean  velocity  is 
everywhere  zero. 

The  case  of  uniform  mean  flow  has  attracted  much  attention  because 
it  is  isotropic.  Naturally,  the  condition  of  isotropy  represents  a  very 
considerable  simplification  of  the  turbulence  problem,  and  has  been  ex¬ 
tensively  studied. 

Unfortunately,  however,  this  degree  of  jimplification  is  excessive 
in  the  present  context,  for  with  isotropic  turbulence  and  zero  shear 
rate  of  the  mean  flow  the  average  turbulent  shear  stresses  are  zero,  and 
there  can  be  no  transfer  of  energy  from  the  mean  flow  to  the  turbulence. 
Hence  the  turbulent  energy  dissipated  into  heat  by  viscous  action  is  not 
replenished,  and  the  turbulent  energy  decays  with  time.  The  turbulence 
therefore  is  not  stationaiy  with  respect  to  time,  as  originally  required, 
(it  becomes  steady  only  in  the  limit  as  zero  amplitude  is  approached. 

*  But  this  limit  is  trivial,  in  that  it  really  amounts  to  the  absence  of 
turbulence.) 

For  non-zero  shear  rate,  however,  the  mean  turbulent  shear  stress  is 
also  non-zero,  and  the  sense  is  such  that  net  energy  is  transferred  from 
the  mean  flow  into  the  turbulence.  Hence  an  equilibrium  is  eventually 
reached  between  the  mean  rate  of  energy  input  and  the  mean  rate  of  viscous 
dissipation.  Consequently,  a  stationary  non-zero  amplitude  of  turbulent 
energy  is  ultimately  established,  as  required. 

In  short,  the  required  flow  field  involves 

(a)  A  steady,  incompressible,  parallel  mean  shear  flow  with  constant 
non-zero  shear  rate,  and  constant  viscosity. 
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(b)  Stationary  and  homogeneous  turbulence  which  is,  however,  aniso- 
t  ropic . 

Now  consider  the  flow  between  two  parallel  walls  of  infinite  extent, 
as  shown  in  Pig.  2.1.  The  walls  move  in  opposite  directions  with  con¬ 
stant  velocities  +  U  as  indicated.  Spacing  between  the  walls  is  2h. 

The  fluid  has  constant  uensity  P  and  constant  kinematic  viscosity  v. 

A  Reynolds  number  may  be  defined  for  this  flow  in  the  form 

,  U  h  x 

"e'C  -ir)  (2J,) 


It  is  a  well  known  principle  of  fluid  mechanics  that  for  any  given 
geometrical  configuration,  there  exists  a  corresponding  critical  Reynolds 
number  which  marks  the  boundary  between  laminar  and  turbulent  flow.  For 
the  present  case,  the  critical  Reynolds  number  is  about  750  (Ref.  l) . 
Below  this  critical  value  the  above  flow  is  laminar  and  the  velocity 
varies  according  to  the  simple  linear  law. 


where 


U  «  n  y 


(2-5) 


0 


a  u 
a  y 


constant  shear  rate 


(2-6) 


At  Reynolds  numbers  above  the  critical,  however,  the  flow  becomes 
unstable  with  respect  to  small  velocity  perturbations,  and  a  stationary 
turbulent  condition  is  soon  established.  In  this  case  the  velocity 
distribution  becomes  strongly  non-linear,  somewhat  as  shown  in  the 
diagram. 

Next,  consider  a  sequence  of  turbulent  flows  of  the  above  type  at 
successively  higher  Reynolds  numbers.  Suppose  this  sequence  is  gener¬ 
ated  by  increasing  both  U  and  h  in  such  a  way  that  the  slope  re- 

mains  finite  in  the  mid -region  near  y  ■  0,  far  from  both  walls.  In  the 
limit  as  Rft  -•  ®,  we  once  again  recover  a  linear  distribution  of  the  mean 
velocity  far  from  the  walls,  but  now  the  flow  is  turbulent.  This  is  the 
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Turbulent 

Velocity 

Profile 


situation  we  wish  to  analyse  (see  Fig.  2.2).  It  represents  perhaps 
the  simplest  conceptual  example  of  stationary,  homogeneous  non-isotroplc 
turbulence.  Analysis  of  this  flow  should  reveal  important  basic  informs- 
tion  concerning  the  irreducible  essence  of  stationary  turbulence. 


3.  DIMENSIONAL  ANALYSIS  OF  THE  BASIC  FLOW 


The  basic  flow  field  described  in  the  previous  section,  and  illustra¬ 
ted  lr»  Fig.  2.2,  is  fully  defined  by  three  characteristic  parameters, 
namely , 

f)  ■  ■  the  constant  shear  rate,  l/sec 

p  -  the  constant  fluid  density,  slugs/ft^ 

g 

v  -  the  constant  kinematic  viscosity,  ft  /sec 


The  present  problem  lies  in  the  realm  of  dynamics.  This  amounts  to 
saying  that  all  significant  quantities  which  occur  in  it  are  reducible 
dimensionally  to  certain  combinations  of  not  more  than  three  fundamental 
dimensions.  These  can  be  force  F,  length  L  and  time  T.  However,  the 
requisite  reference  force,  length  and  time  can  be  expressed,  in  turn, 
by  combinations  of  three  parameters  appropriate  to  the  specific  problem 
under  study.  For  the  present  problem,  the  necessary  and  sufficient 
parameters  are  clearly  0,  p  and  v.  These  three  quantities  suffice  to 
establish  a  system  of  natural  reference  units  in  terms  of  which  all 
other  physical  quantities  can  be  expressed  in  an  appropriate  dimension¬ 
less  form.  Thus,  for  example,  if  X  is  any  relevant  quantity  having 
arbitrary  physical  dimensions,  a  dimensionless  form  X*  may  always  be 
defined  cuch  that 


X* 


p“  nb  vc 


The  appropriate  exponents  a,  b,  c  can  always  be  found  so  as  to  render  X* 
dimensionless. 

This  principle  is  illustrated  by  the  examples  listed  in  Table  3*1. 

The  last  three  entries  in  the  table  show  that  the  basic  dimensional 
reference  quantities,  namely  p,  v  and  Cl,  when  themselves  expressed  in 
dimensionless  terms,  turn  out,  of  course,  to  have  unit  amplitude. 
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TABLE  3.1 


\ 


Summary  of  Dimensionless  Variables 


Quantity 

Symbol 

Dimensions 

Dimensionless 

Form 

Principal  Variables: 

Length  x 

L 

**  - x  yy 

Time 

t 

T 

t*«  ot 

Velocity 

u 

L/T 

u*  =  JL_ 

*Tv n 

Velocity  Time 
Derivative 

3u 

at 

l/t2 

au*  i 

Force  per 

Unit  Mass 

f 

F  L 

x..  f 

Pressure  per 

Unit  Density 

0  -  2 
y  P 

2 

FL  L 

Kinetic  Energy 
per  Unit  Mass 

E 

2 

FL  L 

MV 

**•55- 

Reference  Parameters : 

Density 

P 

M  F?2 

“3‘T 

LJ  L 

P*  ■  ^  ■  1 

Kinematic 

Viscosity 

V 

L2 

T 

„  V 

v*  ■  -  ■  X 

V 

0*  - 


Mean  Shear 
Rate 


n 


1 

T 


0 


1 


cvlcv 


This  fact  16  important.  It  shows  that  all  cases  of  unbounded  parallel 
flow  with  constant  shear  rate  are  reducible,  when  expressed  dimension- 
lessly,  to  a  single  easel  This  represents  a  substantial  simplification 
and  generalization  of  the  problem.  Using  dimensionless  nomenclature,  we 
may  say  that  all  cases  are  equivalent  In  a  certain  sense  to  the  case  of 
a  fluid  of  unit  mass  and  unit  viscosity  undergoing  unit  shear  rate. 

It  can  be  shown  rigorously  that  any  equation  that  is  valid  when  ex¬ 
pressed  in  absolute  dimensional  form  Is  also  valid  when  expressed  in 
corresponding  dimensionless  form.  Hence  we  may  replace  absolute  quanti¬ 
ties  x,  t,  u,  . etc.  by  the  equivalent  dimensionless  forms  denoted 

by  x*,  t*,  u* . etc.  In  this  process  the  quantities  p,  v,  0  are  re¬ 

placed  by  p*  =  1,  v*  =  1,  0*  =  1  and  hence  cease  to  appear  explicitly  in 
the  equations. 

The  above  non-dimensionalizing  procedure  is  followed  in  some  of  the 
mathematical  developments  of  this  report.  Whenever  the  context  is  such 
that  there  i6  no  ambiguity  involved,  it  becomes  permissible  to  drop  the 
symbol*  which  was  introduced  above  to  distinguish  dimensionless  forms 
from  their  dimensional  counterparts.  On  the  other  hand,  in  some  discus¬ 
sions  it  is  clearer  to  retain  the  fully  dimensional  form.  This  variation 
of  usage  should  cause  no  difficulty,  as  the  context  makes  the  intended 
meaning  clear. 
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4.  BLOCK  SIZE,  CELL  SIZE,  AND  MEMORY  SIZE 


In  analysis  of  the  flow  by  finite  difference  methods,  it  is  obviously 
necessary  to  confine  the  numerical  operations  to  a  finite  region.  Exten¬ 
sion  of  the  boundaries  to  infinity,  such  as  is  frequently  convenient 
when  one  is  using  analytical  techniques,  is  not  feasible  when  numerical 
methods  are  employed. 

Fortunately,  an  adequate  approximation  can  be  made  on  the  basis  of 
a  finite  region,  provided  that  the  chosen  region  is  sufficiently  large 
to  constitute  a  representative  sample  of  the  field.  The  term  "large"  is, 
of  course,  relative.  We  must  ask,  "Large  relative  to  what?" 

The  answer  is  that  the  region  should  be  large  relative  to  the  "scale 
of  turbulence."  There  are  various  possible  definitions  of  scale  of 
turbulence.  Perhaps  the  most  useful  definition  in  the  present  context 
is  based  on  the  concept  of  statistical  correlation.  Let  us  therefore 
summarize  the  essential  features  of  this  concept. 

For  this  purpose,  let  Xp  and  ^  denote  the  position  vectors  of  two 
arbitrary  but  known  fixed  points  in  the  flow  field.  The  relative  dis¬ 
placement  vector  between  the  two  points  is  then 


Let  Up(t)  and  u^(t)  be  the  respective  turbulent  velocity  vectors  at  the 
two  points  as  functions  of  time.  Assuming  that  the  turbulence  is  ste 
tionary,  we  may  define  the  following  three  time -average  scalar  product  , 
namely, 

T 

up*uQ  =  |  J  Up(t)*uQ(t)  dt  (4-2) 

0 

_  _  T 

Up*up  =  u|  =  |  Jup(t)*up(t)  dt  (4-3) 

0 


1 6 


(4-U) 


_  __  T 

\\  "  uQ  =  I  I  uQ(t)*uQ(t)  dt 
0 

where  the  total  time  T  is  large. 

Now  the  correlation  coefficient  R  is  defined  as 


)  = 


P  Q 


(^-5) 


Since  the  field  is  also  homogeneous,  this  expression  simplifies 
further  because  in  this  case 

~2  ~2  ~2  ,,  0 

UP  =  UQ  *  u  C1*-®) 

Therefore 

A  A 

UP*UQ 

R(?Pq)  -  -~-S-  (M) 

u 


Furthermore,  in  a  homogeneous  field  R  has  the  same  value  for  every 
possible  choice  of  the  point  xp,  provided  ?  is  held  fixed. 

Now  it  is  easy  to  show  that  if  point  Q  is  chosen  to  coincide  with 
point  P,  we  obtain  perfect  correlation,  that  is, 


R(0)  -  1 


(*-8) 


On  the  other  hand,  we  know  from  experiment  that  if  point  Q  is  remote 
from  P,  the  velocities  are  completely  uncorrelated.  Hence 


R  (•)  -  0  (4-9) 

There  is  much  experimental  information  available  which  shows  that 
the  correlation  coefficient  R  diminishes  rapidly  with  increasing  magni¬ 
tude  of  the  relative  displacement  vector,  r  ■  |  r^l.  In  principle, 
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therefore,  there  exists  some  finite  value  ru.Y  beyond  which  the  magnitude 

i  i 

of  the  correlation  coefficient  |  R  |  is  everywhere  less  than  sane  very 
small  preassigned  quantity  e.  Thus 


R<rMAX>  |5,<<<1 


The  actual  value  of  r. 


MAX 


(u-io) 

depends  on  the  permissible  error  e.  Of 


course,  in  any  practical  case  there  is  always  some  limitation  in  the 
accuracy  with  which  R  itself  can  be  measured  or  computed.  If  we  choose 
€  equal  to  the  uncertainty  in  R  itself,  then  there  exists  some  corres¬ 
ponding  finite  We  can  say  that  beyond  this  distance  the  correla¬ 

tion  R  is  truly  negligible.  Hence  the  smallest  region  vhich  is  suffici¬ 
ently  large  with  respect  to  the  scale  of  turbulence  to  be  acceptable  as 
a  representative  sample  of  the  field  may  be  conceived,  for  example,  as  a 
sphere  of  radius 

Another  way  to  Judge  the  size  of  the  region  is  in  relation  to  the 
wavelengths  of  the  turbulent  velocity  spectral  components.  It  is  known 
that  the  wavelengths  which  play  a  significant  role  in  entraining  energy 
from  the  mean  flow,  in  sustaining  the  kinetic  energy  of  the  turbulence, 
and  in  dissipating  turbulent  energy  into  heat,  are  largely  confined 
within  a  certain  limited  range.  The  flow  region  selected  for  study  must 
be  large  enough  to  include  this  range,  of  course,  but  any  further  in¬ 
crease  in  size  can  be  expected  to  have  negligible  effects  on  the  results. 

Provided  the  region  selected  is  sufficiently  large  along  its  smallest 
dimension,  its  exact  shape  should  make  no  noticeable  difference  in  the 
final  results.  Hence  shape  may  be  chosen  on  the  basis  of  convenience. 

Now  the  shape  which  is  by  far  most  convenient  for  numerical  analysis  is 
a  simple  cubical  block,  say  of  length  L  along  each  edge. 

The  cubical  control  block  must  in  turn  be  subdivided  into  small 
volume  elements  for  purposes  of  numerical  analysis  by  differencing  tech¬ 
niques.  One  again,  the  shape  of  the  volume  elements  is  relatively  unim¬ 
portant,  provided  they  are  reasonably  compact  and  sufficiently  small. 
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Clearly,  it  is  by  far  most  convenient  to  choose  these  as  cubical  cells 
of  length  (l/n)  along  each  edge.  It  is  essential,  however,  that  the 
length  (Ij/n)  be  sufficiently  small  relative  to  the  scale  of  turbulence. 

In  particular,  it  must  be  smaller  than  the  shortest  wavelength  in  the 
significant  wavelength  range  as  described  above. 

It  follows  that,  ideally,  we  should  like  to  use  a  block  length  L 
which  lies  above  some  upper  critical  limit,  say  L  x,  and  simultaneously 
to  use  a  cell  length  (l/N)  which  lies  below  some  lower  critical  limit, 
Lmin*  *n  PrinciPQl>  this  can  he  achieved  by  choice  of  a  sufficiently 
large  value  of  N,  that  is,  II  >  . 

The  difficulty  to  be  faced  is  that  computer  memory  limitations  will 
make  it  impossible  to  use  a  value  of  N  large  enough  to  satisfy  the  fore¬ 
going  requirement.  This  is  because  computer  storage  requirements  increase 

•5 

roughly  in  proportion  to  N  .  Moreover,  computing  tirr.es  increase  even 

I* 

more  sharply,  perhaps  in  proportion  to  N  .  Consequently,  some  type  of 
partial  simplification  of  the  problem  becomes  mandatory. 

One  answer  to  this  difficulty  is  to  split  the  overall  turbulence 
into  several  ranges  of  wavelength,  and  solve  them  in  succession.  The 
general  idea  has  already  been  explained  briefly.  It  is  believed  that 
existing  third -generation  computers  have  now  reached  a  state  of  develop¬ 
ment  such  that  the  problem  can  probably  be  solved  in  just  two  stages, 
namely,  by  division  of  the  turbulence  into  small-scale  and  large-scale 

effects.  We  can  do  this  by  choosing  a  length  L  which  represents  the 

s 

block  size  for  the  small-scale  turbulence,  and  which  is  at  the  same  time 
equal  to  the  cell  size  of  the  large-ecale  turbulence. 

Thus  we  may  write  - 

For  the  small -*cale  turbulence: 


Cell  Size 
Block  Size 


min\ 


L 

pSN 

min 


(4-11) 
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For  the  large-ecale  turbulence: 


Cell  Size 


Block  Size 


L  ) 
8  ) 
Lmax> 


(4-12) 


It  follows  that  for  a  two -stage  analysis  of  this  kind,  the  required 
value  of  N  is  reduced  to 


Another  problem  is  that  the  appropriate  spectral  wavelengths  L 

max 

and  Lm1n  are  initia  y  unknown,  although  existing  experimental  data 
might  possibly  provide  this  information.  However,  we  nov  consider  a 
method  of  numerical  experimentation  on  the  computer  itself  which  can 
serve  to  establish  a  reasonable  estimate  for  the  lower  limit  L-n . 

Recall  that  energy  input  to  the  turbulence  is  mostly  concentrated  at 
the  longer  wavelengths  while  energy  dissipation  is  mainly  in  the  shorter 
wavelength  range.  Hence  a  reduction  of  the  small-scale  block  size  L 

B 

amounts  to  cutting  down  energy  input  in  relation  to  dissipation.  It 

should  result  in  a  reduction  in  the  mean  steady  kinetic  energy  of  the 

turbulence.  In  fact,  it  should  presumably  be  possible  to  find  a  critical 

limit  L  below  which  the  initially  imposed  turbulence  cannot  sustain 
cr 

itself  and  eventually  dies  out.  This  condition  can  be  associated  with 
a  critical  Reynolds  number  marking  the  borcer  between  laminar  and  turbu¬ 
lent  flow.  Thus 


L2  n  2 

-£-r—  .  l*2 

v  cr 


0*-13) 


It  may  be  argued  that  over  distances  smaller  than  Lcr,  significant 
turbulence  effects  cannot  be  sustained.  This  therefore  provides  a 
reasonable  measure  for  the  minimum  necessary  cell  size.  That  is,  for 
the  small-scale  turbulence  we  may  place 


r 


(4-14) 


L 

—  -  L 
N  cr 

Presumably  It  Is  possible  to  find  L  by  successive  trials,  using  only 

cr 

the  basic  turbulence  computer  program  itself.  Let  a  series  of  runs  be 
made  varying  only  the  block  size  L.  For  trial  values  of  L  smaller  than 
Lcr,  the  initial  turbulence  may  be  expected  always  to  die  out.  For 
values  greater  than  Lcf,  the  turbulence  kinetic  energy  may  be  expected 
eventually  to  reach  a  non-zero  mean  stationary  state.  Hence  a  systema¬ 
tic  series  of  trials  should  theoretically  suffice  to  establish  Lcr  to 

within  any  prescribed  margin  of  uncertainty.  When  L  has  been  found  in 

cr 

this  way,  we  may  then  fix  the  required  block  size  for  the  small-scale 
turbulence  as  equal  to 


L  =  N  L  (4-15) 

s  cr  \  si 

Then  the  required  block  size  for  the  large-scale  turbulence  becomes 

L.  =  NL  =  N2L  (4-lc, 

nL  s  cr  '  ' 

If  the  large-scale  block  size  obtained  in  this  way  satisfies  the 
condition 


h  4  t4-^) 

then  all  requirements  relating  to  block  and  cell  size  will  have  been  met. 

Naturally,  the  higher  the  available  value  of  N,  the  more  accurate 
are  the  results  obtained  in  this  way.  Theoretically,  the  exact  result 
corresponds  to  the  unattainable  hypothetical  limit  of  infinite  N.  It  is 
believed,  however,  that  by  use  of  the  two-stage  approach  outlined  above, 
reasonably  accurate  results  can  be  obtained  with  modest  values  of  N, 
well  within  the  memory  capacity  of  modern  computers.  With  computer 
technology  itself  still  advancing  at  a  brisk  rate,  there  is  every  reason 
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to  suppose  that  the  accuracy  and  adequacy  of  the  results  attainable  will 
improve  steadily. 

Of  course,  the  use  of  two  (or  more)  successive  ranges  of  wavelengths, 
although  it  should  enable  us  to  cc  voensate  in  part  for  limitations  of 
computer  capacity,  is  not  without  its  drawbacks.  For  one  thing,  the 
basic  equations  become  more  complex  as  the  quadratic  terms  in  the  momen¬ 
tum  equation  introduce  interaction  effects  among  the  several  wavelength 
regimes.  A  detailed  treatment  of  the  interaction  problem  lies  outside 
the  scope  of  the  present  report.  It  is  believed,  however,  that  this 
aspect  can  be  successfully  treated. 
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5.  EQUIVALENT  REFERENCE  FRAMES 


The  present  analysis  is  concerned  with  homogeneous  turbulence.  This 
means  that  all  statistical  features  of  the  turbulence,  such  as  the  mean 
kinetic  energy,  for  example,  are  uniform  over  the  entire  flow  field. 
Consequently,  it  is  immaterial  Just  where  the  control  block  is  located 
in  the  flow  field.  In  fact,  we  may  subdivide  the  entire  unbounded  flow 
region  into  a  cubical  array  of  blocks,  somewhat  as  shown  in  Fig.  5*1> 
and  select  one  of  these  at  random  as  the  system  to  be  analyzed. 

However,  if  the  various  blocks  are  to  be  considered  equivalent  to 
one  another,  they  must  all  be  related  to  the  mean  flow  in  a  corresponding 
manner.  Consider,  for  example,  the  array  of  blocks  in  Fig.  5*1*  Block 
outlines  are  indicated  by  heavy  lines,  cell  outlines  within  each  block  by 
lighter  lines.  Three  rows  of  blocks  are  shown  in  the  figure.  Each  of 
these  rows  has  a  mid-plane,  which  is  the  horizontal  plane  midway  between 
the  top  and  bottom  boundaries  of  that  row.  The  vertical  spacing  between 
the  successive  mid- planes  is  equal  to  the  block  size  L.  Let  the  middle 
row  of  blocks  be  regarded  as  fixed  in  the  diagram.  Suppose  the  mean 
velocity  of  the  fluid  is 


U  =  i  n  y  (5-1) 

where  y  is  measured  vertically  from  the  mid-plane  of  the  middle  row. 

Suppose  that  the  upper  row  of  blocks  is  moving  to  the  right  with  a 
speed  fiL  ,  and  the  lower  row  is  moving  to  the  left  with  speed  ftL  . 
Clearly  there  is  now  zero  relative  velocity  between  fluid  and  block  at 
the  mid-planes  of  all  three  rows.  An  observer  in  any  one  row  of  blocks, 
moving  with  that  row,  sees  within  his  own  row  a  distribution  of  the 
relative  mean  velocity  which  is  identical  with  that  seen  by  any  other 
observer  in  the  latter's  own  row.  The  velocity  diagrams  illustrate  this 
fact.  This  is  what  is  meant  by  the  statement  that  all  blocks  must  be 
related  to  the  mean  flow  in  a  corresponding  way. 
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(b)  Blocks  and  cells  at  other  times 


Fig.  5.1  Relative  Positions  and  Velocities  of  Equivalent  Reference  Blocks. 

12  3 

(a)  Blocks  and  cells  at  times  t  ■  0,  p,  - 

(b)  Blocks  and  cells  at  other  times. 

2k 
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(c)  Mean  flow  velocities  relative  to  successive  rows  of  blocks 


Fig.  5.1  (cont'd)  (c)  Mean  Flow  Velocities  Relative  to  Successive  Rows 
of  Blocks. 
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Nov  any  one  of  these  blocks  may  be  chosen  at  random  as  the  system 
to  be  analyzed.  From  the  point  of  view  of  a  hypothetical  observer  within 
that  block  and  moving  along  with  it,  the  block  and  its  associated  grid  of 
cells  constitutes  a  stationary  or  Euierian  frame  of  reference.  All  the 
subsequent  analytical  equations  are  to  be  understood  as  being  expressed 
with  respect  to  such  an  Euierian  reference  frame.  It  is  an  important 
point,  however,  that  the  different  rows  of  blocks  are  to  be  considered 
as  lying  in  different,  though  equivalent,  reference  frames.  Since  all 
of  these  reference  frames  are  moving  in  straight  lines  parallel  to  one 
another  and  at  constant  relative  velocities,  they  are  all  inertial 
systems . 

It  is  convenient  to  assume  that  at  some  initial  time  t  =  0,  all  of 
the  blocks  in  the  field  are  arranged  in  a  simpDe  cubical  array,  with 
their  front  and  rear  surfaces  lying  in  common  vertical  planes  through 
all  rows,  in  the  manner  shown  in  Fig.  5*l(a)*  As  time  passes,  the 
separate  rows  become  gradually  more  staggered  relative  to  one  another, 
owing  to  their  relative  motion,  as  shown  in  Fig.  5.1(b).  However,  after 
each  time  interval  of  length  (l/O) ,  the  relative  movement  between  any 
two  successive  rows  equals  exactly  one  additional  block  length,  so  that 
a  simple  cubical  alignment  is  once  again  established.  It  is  seen, 
therefore,  that  the  steady  and  progressive  relative  sliding  of  the 
various  blocks  gives  to  the  structure  of  the  reference  framework  a 
periodically  recurring  character,  with  period  T  •=  l/O.  Expressed  in 
appropriate  dimensionless  time  units,  the  period  becomes  simply  T*  = 
l/O*  =  1. 

There  is,  however,  a  time  span  which  is  a  true  characteristic  of 
the  turbulence,  and  which  should  be  mentioned  at  this  poiit.  It  is 
associated  with  the  concept  of  correlation  over  time,  and  is  entirely 
analogous  to  the  previously  discussed  concept  of  correlation  over 
distance. 

A 

Let  Up(t)  and  up(t+T)  be  the  velocity  vectors  at  a  fixed  point  Xp 
at  times  t  and  (t+T)  respectively.  If  the  turbulence  is  stationary  we 
can  define  the  time  averages 
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up(t)*dp(t+T)  -  |  J^(t).up(t4T)dt 


(5-2) 


u2"  =  Up(t)*up(t)  »  up(t+T)  *up(t+T) 

T  T  T 

|  J  ^p(t)  •  $p(t)dt  =  |  J  Up(t+T)*Sp(t+T)dt 


(5-3) 


0  0 

Now  the  correlation  coefficient  R  is  defined  as 


R(T) 


Ot)-u_(t+T) 


u 


(5-5) 


If  the  turbulence  is  homogeneous,  then  for  a  fixed  T,  the  value  of 

A 

R  is  the  same  for  every  possible  choice  of  the  point  Xp. 

It  is  seen  that  if  we  place  t  ■  0,  we  obtain  perfect  correlation, 
that  is, 

R(0)  -  1  (5-6) 

On  the  other  hand,  we  know  from  experiment  that  at  times  suffici¬ 
ently  far  apart,  the  velocities  are  completely  uncorrelated,  that  is 


R(oo)  b  o 


(5-7) 


Now  it  follows  that  there  exists  some  time  span  t-T  beyond  which 
the  value  of  R  differs  from  zero  by  an  amount  less  than  the  small  inher¬ 
ent  uncertainty  in  R  itself. 

It  follows  that  if  we  wish  to  calculate  any  time -average  property 
of  the  turbulence,  an  integration  over  a  time  interval  of  say  2T  will  be 
ample  to  provide  a  Just  average.  Also,  if  we  wish  to  approximate  the 
turbulence  by  Fourier  methods,  it  is  permissible,  for  example,  to  treat 
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the  phenomenon  as  periodic  in  time,  provided  that  ve  choose  the  basic 
reference  period  as  at  least  equal  to  2T.  Actually,  our  present  method 
of  solution  does  not  assume  this  kind  of  periodicity  in  time;  it  vill 
be  seen,  however,  that  the  assumption  of  periodicity  in  space  proves 
to  be  most  useful. 
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6.  BOUNDARY  AND  INITIAL  CONDITIONS 


Solution  of  the  fluid  turbulence  problem  involves  the  application 
of  the  principles  of  mass  and  momentum  conservation  to  all  the  indivi¬ 
dual  small  cubical  cells  which  comprise  the  large  fluid  block  under 
study.  The  resulting  difference  equations  which  embody  these  physical 
principles  may  be  termed  field  equations  because  they  apply  to  all  the 
small  volume  elements  throughout  the  field.  They  play  a  role  in  the 
numerical  difference  method  analogous  to  that  of  partial  differential 
equations . 

However,  such  a  set  of  field  equations  is  capable  of  establishing 
an  indefinitely  large  number  of  distinct  solutions,  depend- ng  on  the 
boundary  and  initial  conditions  which  are  involved.  For  the  present 
case,  the  boundary  conditions  are  simply  the  velocity  rad  pressure  dis¬ 
tributions  occurring  as  functions  of  time  over  the  six  square  bounding 
surfaces  of  the  block.  The  initial  conditions  are  the  distributions 
existing  over  the  entire  control  volume  at  the  initial  time  t  =  0. 

From  a  purely  mathematical  point  of  view,  the  boundary  and  initial 
conditions  may  be  arbitrarily  prescribed;  the  field  equations  then 
establish  the  corresponding  complete  time -dependent  distributions  over 
the  volume  of  the  block.  Th’xs  the  field  equations  do  not  in  themselves 
fully  determine  the  phenomena,  but  merely  propagate  the  effects  of  the 
boundary  and  initial  conditions  in  the  proper  manner. 

But  for  the  purpose  of  simulating  turbulence,  what  are  the  appropri¬ 
ate  boundary  and  initial  conditions  to  apply?  The  field  equations  them¬ 
selves  do  not  enlighten  us  directly  on  this  crucial  point.  Nevertheless, 
there  are  other  considerations  which  do  lead  to  fairly  satisfactory 
answers.  In  this  connection,  it  is  convenient  to  consider  the  initial 
conditions  and  the  boundary  conditions  separately. 

In  regard  to  the  initial  conditions,  it  is  shown  elsewhere  in  this 
report  that  for  a  simple  cubical  grid  of  N^  points,  there  are  2N^  initial 
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degrees  of  freedom.  This  means  that  2N~  initial  velocity  components  may 
be  prescribed  arbitrarily.  Hence  the  number  and  variety  of  conceivable 
initial  combinations  is  extremely  large.  So  how  shall  a  specific  suit¬ 
able  initial  combination  be  selected  or  prescribed? 

The  key  to  the  answer  lies  in  the  fact  that  fluid  turbulence  i6  an 
instability  phenomenon.  It  is  characteristic  of  instability  processes 
that  initial  perturbations  serve  mainly  as  a  triggering  mechanism;  the 
state  ultimately  attained  is  essentially  independent  of  the  exact  nature 
of  the  triggering  disturbance.  The  reason  is  that  certain  components  of 
the  turbulence  tend  to  grow  and  others  to  decay  until  an  inner  dynamic 
equilibrium  is  attained.  This  ultimate  equilibrium  depends  on  the  intrin 
sic  mechanism  of  turbulence  rather  than  upon  the  exact  form  of  the 
triggering  perturbation. 

It  is  known,  however,  that  the  initial  response  of  a  particular 
flow  field  to  an  imposed  disturbance  may  be  relatively  sensitive  to 
certain  ranges  of  frequencies  or  wavelength  and  comparatively  immune  to 
others,  although  the  ultimate  state  is  independent  of  these  details. 

Hence  it  is  probably  advantageous,  though  not  essential,  to  assume  an 
initial  velocity  perturbation  which  contains  a  wide  range  of  wavelengths 
and  frequencies.  Apart  from  this  rather  mild  requirement,  however,  the 
exact  form  of  the  initial  perturbation  should  be  immaterial  in  the  long 
run. 

Fran  this  it  appears  that  any  random  set  of  initial  velocity  pertur¬ 
bations  is  adequate  for  starting  the  calculation  procedure.  In  any  case, 
the  precise  influence  traceable  to  the  specific  initial  input  used  may 
be  expected  to  fade  out  with  time  in  a  roughly  exponential  manner, 
quickly  leading  to  a  state  of  the  system  which  may  for  all  practical 
purposes  be  regarded  as  stationary  and  representative.  At  least,  this 
is  the  assumption  made  in  the  present  analysis. 

Turning  now  to  the  conditions  along  the  bounding  surfaces,  we  find 
a  more  difficult  problem.  The  reason  is  that  the  influence  of  the 
boundaries,  unlike  that  of  the  initial  state,  does  not  diminish  with 
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time.  Hence  there  is  no  possibility  of  obtaining  a  valid  solution 
unless  the  boundary  conditions  are  defined  in  an  appropriate  way. 

In  this  connection,  one  important  clue  is  that  the  flow  field  is 
homogeneous.  Therefore,  the  exact  placement  of  the  reference  grid  is 
arbitrary  and  immaterial.  Thus  the  overall  statistical  results  obtained 
from  a  certain  system  of  blocks  should  be  the  same  as  those  from  another 
system  of  blocks  displaced  with  respect  to  the  first  by  arbitrary  amounts 
along  all  three  axes.  Of  course,  if  there  is  a  vertical  displacement 
between  these  two  reference  systems,  there  must  also  be  a  corresponding 
horizontal  relative  velocity  between  them .  As  has  already  been  shown, 
this  is  necessary  in  order  that  both  systems  of  reference  blocks  shall 
have  corresponding  velocity  relationships  with  respect  to  the  mean  flow 
passing  through  them. 

Another  way  to  view  the  matter  is  that  the  translation  of  any  plane 
by  any  amount  in  any  direction  should  not  change  any  overall  statistical 
characteristics  pertaining  to  that  plane,  provided  only  that  it  remains 
parallel  to  its  original  orientation  and  that  it  retains  a  corresponding 
velocity  relative  to  the  mean  flow. 

Another  clue  concerning  the  character  of  the  boundary  conditions  is 
that  each  bounding  surface  is  the  common  interface  between  two  adjacent 
blocks.  It  serves  as  a  medium  of  transmission  of  the  influence  of  each 
of  the  blocks  upon  the  other.  Furthermore,  there  is  a  kind  of  symmetry 
in  this  two-way  transmission.  The  reason  is  that  both  blocks  are  suffici¬ 
ently  large  that  all  statistical  characteristics  of  the  turbulence  are 
sensibly  equal  in  the  two  volumes.  It  is  seen  therefore  that  each  bound¬ 
ing  surface  is  the  common  interface  between  tvo  large  volumes  whose  fluid 
contents  exhibit  statistically  equivalent  behavior. 

Another  significant  attribute  of  turbulence  that  should  be  borne  in 
mind  in  formulation  of  appropriate  boundary  conditions  is  the  element  of 
randomness  in  the  phenomena.  The  statistical  concept  of  correlation,  as 
explained  in  an  earlier  section,  is  helpful  in  this  connection.  In 
particular,  we  require  that  the  velocity  correlation  coefficient  R(?pq) 
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shall  be  negligible  If  the  distance  between  the  two  points  involved 

lx 

is  sufficiently  large,  say  equal  to  block  size  L. 

Consider  specifically  two  points  P  and  Q,  point  P  being  somewhere 
on  the  left  bounding  surface  and  point  Q  being  somewhere  on  the  right 
bounding  surface  of  the  block.  The  distance  between  P  and  Q  is  never 
less  than  the  block  size  L  regardless  of  the  precise  locations  of  P  and 
Q  on  the  respective  faces.  Hence  the  velocity  correlation  between  P 
and  Q  must  be  negligible  for  all  possible  pairs  of  locations.  Similar 
conclusions  can  be  drawn  for  another  pair  of  points  P'  and  Q'  on  the 
front  and  rear  bounding  faces  of  the  block,  respectively.  Similarly 
for  points  P"  and  Q"  on  the  upper  and  lower  faces. 

At  first  glance,  the  above  restriction  seems  to  rule  out  the  idea 
of  treating  the  spacewise  velocity  components  as  periodic  functions  of 
wavelength  L,  equal  to  block  size.  For  if  such  periodicity  were  assumed 
then  pairs  of  points  could  always  be  found  on  opposite  faces  for  which 
the  correlation  coefficient,  instead  of  vanishing,  would  equal  unity! 

What  seems  to  be  demanded,  instead,  is  a  certain  statistical  ran¬ 
domness  in  the  boundary  conditions.  It  is  well  understood,  of  course, 
that  turbulence  is  a  partly  stochastic  process.  Since  the  field  equa¬ 
tions  themselves  are  strictly  deterministic,  the  stochastic  aspects  can 
only  enter  via  the  initial  and  boundary  conditions. 

It  is  undoubtedly  possible  to  devise  some  calculation  procedure  of 
tte  Monte  Carlo  type  to  simulate  the  stochastic  element  in  the  boundary 
conditions,  while  still  preserving  the  various  types  of  statistical  equi 
valence  previously  discussed.  The  advantage  of  such  a  procedure  would 
be  that  all  the  foregoing  requirements  could  be  met,  including  the  re¬ 
quirement  for  high  correlation  at  small  distances  and  zero  correlation 
at  large  distances.  Hence  thiB  possibility  is  worth  more  detailed 
study. 

It  has  been  concluded,  however,  that  a  very  slight  relaxation  of 
the  above  requirements  leads  to  a  great  simplification  of  the  problem, 
without  any  serious  disadvantages.  Also,  It  effectively  minimizes  the 
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stochastic  aspects,  and  puts  the  calculations  on  an  essentially  determin¬ 
istic  basis. 

The  simplification  consists  of  selecting  for  the  initial  perturbation 
a  pattern  which,  although  it  may  be  somewhat  random  in  other  respects, 
is  "blockwise  periodic."  By  this  term  we  mean  that  at  the  initial 
instant  of  time,  the  distributions  of  the  turbulent  velocity  and  pressure 
perturbations  are  identical  in  all  blocks. 

Now  it  follows  rigorously  from  the  equations  of  motion  that  if  the 
turbulence  is  blockwise  periodic  at  the  initial  instant,  it  remains 
blockwise  periodic  thereafter. 

To  define  this  concept  mathematically,  let  us  take  the  typical 
scalar  velocity  component  u(x,y,z,t),  for  example.  The  function  u  will 
be  said  to  be  blockwise  periodic  if,  and  only  if,  for  all  values  of 
x,y,z  and  t  we  have 


u(x*,  y',  z',  t)  *  u(x,y,z,t) 


where  x1  «=  x  +  pL  +  qLfit  ) 

y*  *  y  +  qL  ) 

z'  a  z  +  rL  ) 


and  where  p  ) 

q  )  *  0,  +  1,  +  2,  +  3,  ••••  +  00 

r  )  -  -  - 


(6-1) 


(6-2) 


(6-3) 


The  term  of  qLOt  in  the  above  expression  accounts  for  the  difference 
in  mean  flow  velocity  in  two  horizontal  planes  separated  by  the  vertical 
distance  qL. 

This  assumption  of  periodicity  among  the  blocks,  is,  of  course, 
stronger  than  the  requirement  for  merely  statistical  equivalence.  The 
distributions  in  the  various  blocks  could  be  statistically  equivalent 
without  necessarily  being  identical  at  every  corresponding  point  of 
space  and  time.  On  the  other  hand,  if  they  are  treated  as  stiictly 
periodic  in  this  special  way,  then  they  satisfy  statistical  equivalence 
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identically.  That  is  to  say  all  volumetric  integrals,  averages,  vari¬ 
ances  etc.  will  be  exactly  equal  for  all  blocks. 

It  turns  out  that  this  assumption  of  blockwise  periodicity  supplies 
exactly  the  information  needed  to  establish  the  unknown  boundary  condi¬ 
tions  in  a  completely  determinate  manner,  without  recourse  to  auxiliary 
statistical  hypotheses  or  ad  hoc  arguments.  In  other  words,  the  hypo¬ 
thesis  is  very  attractive  because  it  succeeds  so  well! 

The  method  is  also  efficient  numerically  in  the  sense  that  the 
stored  numerical  information  now  is  representative  not  Just  of  one  block 
but  of  all  blocks,  and  hence  of  the  entire  infinite  field. 

One  requirement  that  the  assumption  of  blockwise  periodicity  does 
not  fully  meet  is  that  velocity  correlations  should  vanish  at  large  dis¬ 
tances.  Instead,  it  is  found  that  the  correlation  coefficient  itself 
becomes  periodic  in  character.  However,  there  is  no  great  harm  in  this, 
provided  that  the  basic  wavelength  is  sufficiently  large  that  the  recurring 
correlation  peaks  are  well  separated.  In  fact  if  the  correlation  is 
negligible  at  the  midpoints  between  successive  peaks,  each  peak  may  be 
regarded  as  essentially  isolated,  and  the  spurious  correlation  at  the 
full  wavelength  may  be  safely  disregarded.  Figure  6.1  may  help  to 
clarify  these  ideas. 
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7.  NUMERICAL  REPRESENTATION  OP  AN  ARBITRARY 
SCALAR  FUNCTION  OP  POSITION 


In  this  analysis  ve  deal  with  scalar  functions  3uch  as  the  velocity 
components  u,  v,  w,  and  the  pressure  function  cp.  Typically  these  quanti¬ 
ties  are  continuous  functions  of  the  spatial  coordinates  x,  y,  z  and  of 
time  t.  Symbolically  we  may  write,  for  example, 

u  -  u(x,y,*,t)  (7-1) 

Furthermore,  the  various  partial  derivatives  of  all  orders  exist  and 
a?:«  continuous  functions  of  space  and  time. 

At  any  specified  Instant  of  time,  however,  the  quantities  of  inter¬ 
est  are  functions  of  the  space  coordinates  only.  Thus 

u  -  u(x,y,z)  (7-2) 

Unfortunately,  the  functions  of  interest  possess  very  intricate 
spatial  distributions  which  fluctuate  erratically  with  the  passage  of 
time.  In  principle,  these  functions  may  theoretically  be  approximated 
by  Fourier  series  (or  Fourier  integrals).  In  practice, however,  the 
number  of  terms  required  in  the  series  is  of  the  same  order  as  the 
number  of  grid  points  in  the  block.  It  is  easy  to  see  that  round-off 
errors,  even  if  negligible  on  individual  terns  of  the  series,  will 
accumulate  when  summed  over  the  very  large  number  of  terms  required  in 
the  whole  series.  This  leads  to  excessive  errors  which  vitiate  the 
calculation.  The  plain  fact  is,  therefore,  that  we  do  not  have  available 
any  suitable  analytical  expressions  for  describing  such  intricate  dis¬ 
tributions  in  a  way  that  is  adequate  for  purposes  of  numerical  analysis. 

The  problem,  therefore,  is  to  devise  a  numerical  scheme  for  approxi¬ 
mating  the  required  functions  to  the  requisite  degree  of  accuracy.  The 
scheme  must  also  be  flexible  and  efficient. 
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For  this  purpose  we  define  a  suitable  gridwork  of  reference  points 
uniformly  distributed  throughout  the  space  in  question.  The  values  of 
the  function  are  then  specified  for  all  points  of  the  grid.  In  addi¬ 
tion,  a  suitable  three-dimensional  interpolation  rule  is  needed  for 
estimating  the  value  of  the  function  at  an  arbitrary  point  lying  any¬ 
where  in  the  space  between  the  specified  grid  points.  Since  the  func¬ 
tion  being  approximated  is  single-valued  everywhere,  it  is  desirable 
that  the  Interpolation  rule  adopted  be  definite  and  unambiguous.  In 
other  words,  once  the  values  of  the  function  are  specified  at  the 
reference  grid  points,  the  function  should  be  uniquely  determined  over 
the  entire  space. 

In  the  course  of  this  research,  several  different  schemes  of 
differencing  and  interpolation  have  been  conceived.  Three  of  these 
methods  are  described  in  the  succeeding  sections  of  this  report.  They 
are  discussed  in  order  of  increasing  complication.  However,  the  third 
(Method  0)  was  actually  developed  first  and  used  in  the  work  to  be  dis¬ 
cussed.  It  was  only  later  that  it  was  realized  that  the  simple  Method 
A  could  be  made  to  yield  precise  results. 
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8.  DIFFERENCING  AND  INTERPOLATION:  METHOD  A 


Consider  a  cubical  region  of  length  L  on  a  side.  Let  each  side  be 

1 3 

subdivided  into  N  equal  intervals,  thus  defining  a  set  of  NJ  cubical 
cells,  each  cell  of  length  (l/n),  as  illustrated  in  Fig.  8.1.  Let 
these  cells  represent  the  control  volumes  for  the  finite-element 
analysis. 

Each  of  the  above  control  volumes  has  a  definite  centroidal  point. 
The  centroids  of  the  control  volumes  constitute  a  cubical  array  of 

3 

grid  points.  There  are,  of  course,  NJ  grid  points  in  all. 

For  purposes  of  the  present  discussion,  let  us  choose  an  origin  of 
coordinates  as  shown  in  Fig.  8.1.  Then  the  coordinates  of  the  centroi¬ 
dal  points  may  be  expressed  in  the  form 


x 


i 


y 


J 


(1  -1) 

L 

N 

i  =  1,2,3,. 

...N 

(8-1) 

(J-i> 

L  L 

N  '  2 

J  =  1,2,3, • 

...N 

(8-2) 

(k  - 1) 

L 

N 

k  =  1,2,3,. 

...N 

(8-3) 

We  wish  to  describe,  to  within  a  suitable  degree  of  approximation, 
the  distribution,  over  the  volume  of  the  block,  of  an  arbitrary  scalar 
function  of  position,  say  u(x,y,z). 

The  first  part  of  the  description  consists  of  specifying  the  values 

3 

of  u  *  u(i,j,k)  at  each  of  the  foregoing  grid  points.  So  far  as 

3 

the  immediate  discussion  is  concerned,  these  N  values  may  be  regarded 
as  arbitrary. 

The  second  pert  of  the  description  consists  of  a  suitable  interpola¬ 
tion  rule  for  determining  the  value  of  the  function  at  any  point  other 
than  a  centroidal  grid  point.  For  this  purpose  we  define  a  second  set 
of  cubical  cells  which  we  may  designate  as  interpolation  spaces  or  cells. 
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Fig.  8.1  Basic  Grid  for  Differencing  Method  A. 
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An  interpolation  space  is  also  a  small  cubical  cell  of  length  L/N  on  a 
side.  However,  it  is  displaced  with  respect  to  the  control  volumes. 

The  previously  defined  grid  points  which  are  at  the  centroids  of  the 
control  volumes  lie  at  the  comers  of  the  interpolation  spaces.  It 
follows  that  the  control  volumes  and  the  interpolation  cells  are  mutu¬ 
ally  staggered  relative  to  one  another.  Thus  each  octant  of  a  contro] 
volume  is  in  a  different  interpolation  space.  Conversely,  each  octant 
of  an  interpolation  space  lies  in  a  different  control  volume. 

Next,  we  adopt  the  rule  that  within  each  interpolation  space,  the 
variation  of  the  function  u  shall  be  linear  along  every  line  parallel 
to  one  of  the  coordinate  axes.  This  rule  can  be  expressed  in  algebraic 
form. 

For  this  purpose,  it  is  convenient  to  introduce  the  following 


auxiliary  notation .  Let 

E  («i) 

<1 

(8-4) 

Tl  = 

E  ("j) 

0  S  T]  <  1 

(8-5) 

C  « 

I  U) 

0  s  C  *  1 

(8-6) 

Thus  the  dimensionless  coordinates  ? ,  T|,  C  i&&y  each  vary  over  the  range 
from  zero  to  unity  within  any  interpolation  space. 

Now  the  variation  of  u  within  any  particular  interpolation  space 
may  be  expressed  by  the  polynomial 

u  -  aq  +  Ax  ?  +  a^  +  a3c  +  A^ttC  +  a5CS  +  +  A^tiC  (8-7) 

The  A's  are  initially  unknown  constants  whose  values  must  be  determined 
from  the  known  values  of  u  at  the  eight  corners  of  the  particular  inter¬ 
polation  space  in  question.  At  these  eight  corners,  each  of  the  coordi¬ 
nates,  ?,  ti,  £,  is  either  zero  or  unity.  Hence  the  foregoing  equation 


ho 


can  be  written  eight  times,  with  the  appropriate  values  of  T|,  C >  and 
u  substituted  each  time.  The  resulting  set  of  eight  simultaneous  equa¬ 
tions  Just  suffices  to  determine  the  eight  initially  unknown  constants. 
Because  of  the  symmetry  and  the  many  zeros,  the  solution  turns  out  to 
be  quite  simple. 

Note  that  Eq.  (8-7)  does,  in  fact,  satisfy  the  requirement  that  u 
shall  vary  linearly  along  any  line  parallel  to  any  one  of  the  coordinate 
axes.  This  is  because  the  coordinates  §,  r|,  £  occur  individually  only 
to  the  zeroth  or  first  powers.  Furthermore,  this  equation  contains  all 
eight  of  the  possible  combinations  of  these  factors,  and  is  therefore 
the  most  general  possible  expression  capable  of  satisfying  this  three- 
dimensional  linearity  rule. 

Note  that  the  interpolation  spaces  along  the  six  surfaces  of  the 
block  extend  into  the  adjacent  blocks.  However,  the  variation  within 
these  spaces  can  be  handled  in  much  the  same  way  as  for  the  interior 
spaces.  In  this  connection,  it  is  convenient  to  take  advantage  of  the 
aseumed  spacewise  periodicity  of  the  function  u  in  the  x  and  z  direc¬ 
tions.  Thus 

u(i  +  N,j,k)  -  u(i,J,k)  (8-8) 

u(i, j,k  +  N)  =  u(i,j,k)  (8-9) 

There  is  also  a  periodicity  in  the  y  direction,  but  it  is  of  a 
slightly  more  complicated  kind.  This  matter  is  discussed  elsewhere  in 
this  re]>ort.  It  suffices  here  to  say  that  an  appropriate  periodicity 
rule  is  available  for  the  y  direction  which  makes  the  variation  of  u 
perfectly  determinate  also  for  the  interpolation  spaces  along  the  top 
end  bottom  surfaces  of  the  block. 


The  foregoing  diecussion  shows  that  once  the  values  of  u(i,J,k)  are 

•3 

arbitrarily  specified  at  the  NJ  grid  points,  the  above  interpolation 
rule  defines  the  function  u(x,y,z)  uniquely  over  the  entire  volume  of 
the  block. 

This  means,  of  course,  that  the  distribution  of  u(x,y,z)  is  also 
uniquely  determined  over  each  of  the  six  square  surfaces  which  enclose 
each  control  volume.  This  fact  is  important,  because  the  application 
of  the  fundamental  mass  and  momentum  conservation  laws  involves  the 
evaluation  of  certain  integrals  over  these  bounding  surfaces. 

In  addition,  the  determination  of  certain  overall  flow  characteris¬ 
tics,  such  as  the  mean  kinetic  energy  of  turbulence,  for  example,  in¬ 
volves  the  evaluation  of  corresponding  volume  integrals.  Consequently, 
a  unique  rule  for  the  variation  of  any  function  over  the  entire  volume 
is  necessary. 

2  2  2 

The  six  quadratic  velocity  products  u  ,  v  ,  w  ,  vw,  wu,  uv  play  an 
important  role  in  the  analysis.  In  dealings  with  these  products,  the 
usual  quasi-iinear  interpolation  rule  is  taken  to  apply  to  the  indivi¬ 
dual  factors,  rather  than  simply  to  the  product  itself.  Hence  the 
variation  of  the  product  between  grid  points  is  quadratic  along  lines 
parallel  to  the  coordina-e  axes,  rather  than  linear.  While  this  compli¬ 
cates  the  calculations  somewhat,  it  yields  a  higher  level  of  consistency 
and  accuracy  than  could  otherwise  be  attained.  In  view  of  the  basic 
importance  of  these  products,  which  are  associated  with  the  Reynolds 
stresses,  the  added  accuracy  appears  to  be  well  worth  the  extra  complies^ 
tions. 

Since  the  entire  spacewise  distribution  of  any  function  is  deter- 

3 

mined  when  the  values  sre  known  at  the  N  grid  points,  it  follows  that 
all  pertinent  surface  and  volume  integrals  must  be  ultimately  expres¬ 
sible  directly  in  terms  of  these  grid  values.  When  these  relationships 
are  expressed  in  explicit  form,  they  constitute  the  "differencing 
formulas"  appropriate  to  the  task  at  hand. 
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9.  DIFFERENCING  AND  INTERPOLATION:  METHOD  B 

In  this  scheme,  two  distinct  families  of  grid  points  are  used,  as 
shown  in  Fig.  9*1*  It  is  convenient  to  designate  these  families  by- 
means  of  superscripts  and  ,  respectively. 

Points  of  family  I  are  defined  by  the  coordinates 


I 

xi = 

(i  -  §) 

L 

N 

i  =  1,2,3, . 

.N 

(9-1)  I 

I 

yJ = 

(J  -  §) 

L  L 

N  ‘  2 

J  =  1,2,3, . 

.N 

(9-2) 

I 

\  = 

(k  -  §) 

L 

N 

k  =  1,2,3, . 

.N 

(9-3) 

Points  of  family  II  j 

are  defined 

by  the  coordinates 

x11  « 
i 

k|) 

i  =  0,1,2, . 

.N 

(9-4) 

II 

V  .  = 

L 

1  b  0.1.2 . 

.N 

i  •  1 

(9-5) 

2 

II 

z  — 

*<§> 

k  -  0.1.2 . . 

.N 

(9-6) 

1 

k 

Note  that 

the  two  grids  are 

mutually  staggered 

relative 

1 

to  each 

other  by  the  distance  of  half  a  cell  width,  L/2N,  with  respect  to  each 
of  the  coordinate  planes. 

Corresponding  to  the  two  families  cff  interlaced  grid  points,  there 
are  two  families  of  interlaced  control  volumes.  All  control  volumes  are 
cubical  cells  of  length  (L/N)  on  a  side.  The  grid  points  of  family  I 
are  the  centroids  of  the  control  volumes  of  the  first  family.  Likewise, 
the  grid  points  of  family  II  are  the  centroids  of  the  control  volumes 
of  the  second  family.  Furthermore,  the  grid  points  of  family  II  lie  at 
the  corners  of  the  control  volumes  of  family  I.  Similarly,  the  grid 
points  of  family  I  lie  at  the  corn  -s  of  the  control  volumes  of  family  II. 
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These  various  relations  sre  inherently  three-dimensional  and  are 
not  easy  to  portray  adequately  in  two-dimensional  diagrams.  Neverthe¬ 
less,  Pigs.  9-1  and  9.2  may  be  of  some  help  in  visualizing  the  essential 
facts. 

Note  that  the  overall  volume  of  the  cubical  block  under  study  is 
represented  twice  in  this  scheaet  once  by  the  control  volumes  of  family 
I,  and  once  by  the  control  volumes  of  family  II.  There  is  a  mutual 
overlapping  of  the  individual  control  volumes.  Thus  the  individual 
octants  of  a  single  cell  of  family  I  lie  within  eight  separate  cells  of 
family  II.  Conversely,  the  individual  octants  of  a  cell  of  family  II 
lie  within  eight  separate  cells  of  family  I. 

Consider  any  one  such  octant  as  shown,  for  example,  in  Fig.  9*2. 

It  is  itself  a  cube  of  length  (L/2N)  on  a  side.  One  corner  of  this  cube 
is  always  a  grid  point  of  family  I.  The  diagonally  opposite  corner  is 
always  a  grid  point  of  family  II.  We  define  each  distinct  octant  as  a 
distinct  interpolation  space.  There  are  therefore  8n  distinct  inter¬ 
polation  spaces  within  the  block. 

We  note  for  comparison  that  Method  A  involves  only  NJ  distinct  inter¬ 
polation  spaces.  This  does  not  imply,  however,  that  Method  B  involves 
an  eightfold  improvement  in  resolution.  The  reason  is  that  the  inter¬ 
polation  spaces  of  Method  A  are  "complete"  in  the  sense  that  all  eight 
of  their  corners  are  principal  grid  points.  The  interpolation  spaces  of 
Method  B  are  initially  "incomplete"  in  the  sense  that  only  two  of  the 
eight  corners  are  principal  grid  points.  To  "complete"  these  spaces, 
we  must  define  the  values  of  the  unknown  function  at  the  other  six  cor¬ 
ners  of  each  space. 

Fortunately,  each  of  these  six  corners  is  itself  a  midpoint  between 
two  principal  grid  points.  Hence  using  linear  interpolation  consistently 
throughout,  we  may  assign  a  value  at  each  of  these  six  points.  The 
assigned  value  is  always  the  simple  arithmetic  mean  of  the  known  values 
at  the  two  corresponding  principal  grid  points  involved.  Thus,  for 
example,  in  Fig.  9*2  the  value  at  point  QII  is  the  mean  of  the  values  at 

^5 


points  PII  and  All.  Similarly,  the  value  at  point  Q'l  i6  the  mean  of 
the  values  at  points  P'l  and  A'l.  Note  that  three  of  the  interpolated 
values  are  derived  from  values  at  grid  points  of  family  I,  the  other 
three  from  grid  points  of  family  II.  Thus  both  families  are  represented 
exactly  equally  and  symmetrically  within  each  interpolation  space. 

With  the  eight  oorner  values  defined,  we  can  once  more  write  the 
variation  over  the  volume  and/or  surfaces  of  the  interpolation  space  in 
the  form 


u  "  A0+A15+A2n+A3c+A4ric+i\5C?+A6?n+A7?TTC  (9-7) 

where,  for  example, 


,  x  2N 

(x^)  r 

0  *  %  *  1 

(9-8) 

(y-yj)  r 

0  s  n  <  1 

(9-9) 

,  v  2N 

0  <  c  *  1 

(9-10) 

and  where  the  point  at  the  "inner"  comer  x^,  y^,  z^  may  be  of  either 
family,  depending  on  the  particular  interpolation  space  being  considered. 

It  is  apparent  that  the  foregoing  definitions  suffice  to  determine 
uniquely  the  complete  spacewise  distribution  of  any  scalar  function 
u(x,y,z)  in  terms  of  the  assigned  values  at  the  2NJ  principal  grid 
points.  Consequently  all  surface  and  volume  integrals  over  the  control 
cells  can  be  uniquely  expressed  in  terns  of  the  appropriate  grid  point 
values.  Hence  the  method  leads  to  specific  "differencing  formulas"  for 
such  quantities  as  partial  derivatives,  gradients,  divergences,  etc. 
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10.  DIFFERENCING  AND  INTERPOLATION:  METHOD  C 


The  differencing  scheme  described  under  this  heading  is  merely  an 
earlier  and  somewhat  less  fully  developed  version  of  the  staggered 
scheme  described  as  Method  B.  It  ie  important  however,  because  the 
currently  available  computer  program,  TURBOCODE  MARK  V,  is  based  on 
this  system. 

Method  C  uses  the  two  families  of  grid  points  and  control  volumes 
described  previously.  The  interpolation  spaces  and  interpolation  rules 
are  defined  differently,  however.  In  fact  the  interpolation  spaces  in 
this  scheme  are  taken  as  identical  with  the  respective  control  volumes 
themselves;  there  are  therefore  two  distinct  and  overlapping  sets  of 
interpolation  spaces. 

Thus  within  the  volume  and  along  the  surfaces  of  any  cell  of 
family  I,  the  variation  of  any  scalar  function  u(x,y,z)  is  assumed  to 
follow  the  rule 

u1  =  A*  +  Ai?  +  *2*  +  A3C  +  Aifa  +  A5C?  +  A6?T1 


+  A^t<  + 

A^(1-??)(1-T12)(1-C2) 

(10-1) 

?  -  (x-xj)  f 

-i  <  5  <  +i 

(10-2) 

"  ■  <*-p>  r 

-l  s  n  *  +i 

(10-3) 

oilk-j 

M 

1 

M 

II 

vj- 

■H  (  i  +1 

(10-4) 

Similarly  for  any  cell  of  family  II,  the  corresponding  rule  is 
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where 


u11  =  aJZ+  +  A^T,  +  A^c  +  aJ1*  +  A°c?  +  AgX?n 

+  A^SnC  +  AgI(i-^)(i-ii2)(i-c2)  (10-5) 

?  =  (x-xj1)  ~  (10-6) 

ti  =  (y-y*1)  r  -i  *  r\  *  +i  (10-7) 

C  =  (z-z*1)  f*  -1  <  C  5  ■»!  (10-8) 

There  are  a  number  of  critical  observations  to  be  made  about  thiB 
scheme . 

Note  first  that  any  given  arbitrary  point  in  the  field  is  simul¬ 
taneously  in  two  different  interpolation  spaces.  Two  distinct  inter¬ 
polation  formulas  apply  which  do  not,  in  general,  agree  exactly.  Such 
an  ambiguity  is  not  desirable,  although  it  may  be  argued  that  the  dis¬ 
crepancies  involved  are  small  quantities  of  higher  order,  and  hence 
relatively  unimportant.  However,  the  existence  of  two  distinct  rules 
makes  it  necessary  to  compute  certain  volume  integrals  twice,  thus  in¬ 
creasing  the  corresponding  computation  times  in  an  undesirable  way. 

Note  secondly  that  the  interpolation  polynomial  contains  not  eight 
but  nine  constants.  These  are  determined  by  matching  the  known  values 
at  the  eight  corners  of  the  cell  and  the  ninth  known  value  at  the 
centroid . 

It  may  be  seen  that  the  last  term  of  the  interpolation  expression, 
unlike  the  others,  is  quadratic  rather  than  linear  in  the  coordinates. 
The  presence  of  this  term  markedly  complicates  the  computation  of  volume 

integrals  for  the  cell.  This  is  especially  true  in  connection  with  the 

2 

volume  integration  of  quadratic  terms  like  u  ,  uv,  etc.  The  linear 
texms  among  themselves  have  certain  properties  of  symmetry  and  ortho¬ 
gonality  which  tend  to  simplify  the  resulting  expressions.  The 

h9 


characteristics  of  the  quadratic  term  are  more  complex,  and  introduce 
various  additional  terms  into  the  final  results. 

Fortunately,  the  quadratic  term  has  been  so  devised  that  it  vanishes 
at  all  six  of  the  bounding  surfaces  of  the  cells.  Consequently  it  does 
not  enter  into  the  evaluation  of  any  of  the  surface  integrals. 

As  compared  with  Method  B,  it  would  seem  that  Method  C  does  not 
intermesh  the  data  of  the  two  families  as  closely  and  completely  as 
might  be  desirable.  Thus  in  Method  B,  four  of  the  data  points  in  each 
interpolation  space  are  of  cne  family  and  four  are  of  the  other  family. 

In  Method  C,  on  the  other  hand,  eight  of  the  data  points  are  always 
of  one  family,  and  only  one  is  of  the  opposite  family;  furthermore,  the 
influence  of  even  this  one  point  of  the  opposite  family  vanishes  along 
the  bounding  surfaces  of  the  cell. 

Fortunately,  Method  C  doe6  possess  some  compensating  computational 
advantages.  For  one  thing,  the  distribution  of  any  function  u(x,y,z) 
over  any  bounding  surface  is  expressed  solely  in  terms  of  the  grid 
values  at  the  four  corners  of  that  surface.  While  this  assumption  does 
not  fully  utilize  all  of  the  information  actually  available,  it  does, 
of  course,  shorten  and  simplify  the  evaluation  of  all  surface  integrals. 
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11.  COMPARISON  OF  DIFFERENCING  METHODS 


■ 


■ 


In  the  course  of  the  present  research  there  has  been  a  steady- 
improvement  in  the  conceptual  basis  for  the  numerical  differencing 
techniques  employed.  At  first  the  equations  to  be  solved  were  regarded 
simply  as  partial  differential  equations.  The  numerical  techniques  at 
this  stage  consisted  mainly  of  replacement  of  infinitesimals  by  small  but 
finite  differences,  with  a  largely  intuitive  approach. 

It  was  soon  realized,  however,  that  much  better  accuracy  and 
computational  stability  could  be  achieved  for  a  given  mesh  size  by 
adoption  of  a  finite -element  approach.  An  essential  point  is  that  the 
small  but  finite  control  volumes  are  not  treated  as  mere  infinitesimals. 
Accuracy  is  greatly  improved  by  estimation  of  mass  and  momentum  fluxes 
not  merely  from  approximate  local  point  values  but  as  true  surface 
integrals,  in  accordance  with  the  previous  carefully  established  inter¬ 
polation  rules. 

It  was  at  about  this  stage  that  the  interpolation  method  C  was 
formulated  and  successfully  applied  in  the  program  TURBOCODE  MARK  V. 

When  first  formulated,  it  represented  a  considerable  advance  in  numerical 
technique. 

Further  theoretical  progress  since  that  time  makes  it  clear  that 
certain  revisions  of  that  scheme  are  desirable.  These  are  described 
under  Method  B. 

However,  the  same  concepts  that  lead  to  the  above  revision,  when 
followed  to  their  logical  conclusions,  lead  even  further  and  so  to 
Method  A. 

The  main  advantage  of  Method  A  lies  in  its  relative  simplicity. 

It  avoids  the  complexities  of  dealing  with  two  distinct  families  of 
points.  It  avoids  the  duplication  of  covering  the  same  overall  volume 
with  two  distinct  families  of  control  cells.  The  application  of  the 
proper  boundary  conditions  should  be  correspondingly  simplified.  It 
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seecis  that  the  programming  can  be  simplified  and  computation  times 
reduced  by  adoption  of  this  scheme. 

It  is  intended  that  Method  A  vill  ultimately  be  programmed  for  the 
computer.  Meanwhile ,  however,  the  present  program  based  on  Method  C 
represents  investment  of  a  considerable  research  effort.  It  is  believed 
that  this  program,  which  has  only  recently  been  developed  into  a  prac¬ 
ticable  form,  is  capable  of  generating  much  useful  basic  information. 
Hence  the  immediate  effort  will  be  to  obtain  this  information  with  the 
existing  program,  and  to  study  and  evaluate  it,  before  proceeding  with 
any  major  revision  of  the  computer  program. 


12.  INITIAL  NUMBER  OF  DEGREES  OF  FREEDOM 


In  this  discussion,  let  G  represent  the  number  of  grid  points  in 

7 

one  block.  For  differencing  method  A,  G  equals  N  ;  for  differencing 
methods  B  and  C,  G  equals  2N  . 

At  the  initial  instant  of  time,  the  velocity  vector  distribution 
over  the  cubical  block  is  fixed  when  the  distributions  of  the  three 
individual  scalar  components.,  u,  v,  w,  are  fixed.  But  if  we  are  to 
establish  the  initial  distribution  of  u,  values  must  be  specified  at 
all  G  grid  points.  Similarly  for  v  and  w.  Hence  3^  values  are  re¬ 
quired  to  describe  the  initial  velocity  field. 

However,  if  the  fluid  is  incompressible,  not  all  of  these  30 
values  may  be  specified  independently.  In  fact,  we  must  satisfy  a  mass 
conservation  condition  or  continuity  law  for  each  one  of  the  G  control 
volumes  centered  at  each  grid  point.  Hence  there  are  G  constraint 
conditions  to  be  satisfied  among  the  30  velocity  components.  Conse¬ 
quently,  only  2G  of  these  quantities  may  be  independently  prescribed. 

We  say  that  there  are  2G  initial  degrees  of  freedom. 

For  example,  if  the  u  and  v  components  are  arbitrarily  specified 
at  all  grid  points,  the  w  components  are  then  fixed  by  the  foregoing 
continuity  conditions.  If  the  distributions  were  not  blockwise  periodic, 
continuity  constraints  would  determine  w  only  to  within  an  arbitrary 
additive  function  of  x  and  y  alone.  However,  the  assumption  of  period¬ 
icity  eliminates  this  arbitrary  element  and  renders  the  solution  for  w 
unique . 

The  initial  values  of  the  pressures  do  not  count  as  initial  degrees 
of  freedom.  The  reason  is  that  once  the  velocities  are  specified,  the 
pressures  are  also  determined  thereby.  The  pressures  cannot  be  assigned 
independently. 
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It  f ollow8  that  for  differencing  method  A,  there  are  2N  initial 
degrees  of  freedom;  for  differencing  methods  B  and  C,  there  are  UN-' 
initial  degrees  of  freedom. 


However,  the  apparent  arbitrarineaa  of  the  2G  Initial  degreea  of 
freedom  require  a  qualification,  as  it  car  be  somewhat  misleading.  This 
enormous  degree  of  latitude  is  permissible  only  because,  at  the  initial 
instant,  the  turbulence  has  not  yet  attained  its  ultimate  stationary 
and  homogeneous  statistical  stiucture.  Once  this  state  is  reached,  it 
will  impose  certain  corresponding  relations  of  statistical  correlation, 
and  thereby  reduce  the  number  of  arbitrary  or  random  degrees  of  freedom 
remaining.  Nevertheless,  turbulence  necessarily  involves  some  irre¬ 
ducible  core  of  randomness,  and  hence  there  will  always  be  some  large 
residual  number  of  arbitrary  or  random  degrees  of  freedom.  If  we  knew 
what  the  essential  constraints  were,  it  might  be  useful  to  impose  then 
at  the  outset,  much  as  we  have  done  with  regard  to  the  continuity  con¬ 
dition.  But  we  do  not  know  them,  at  least  not  initially,  so  the  matter 
is  academic.  It  is  rather  lucky  that  the  calculation  can  proceed  from 
such  arbitrary  initial  conditions.  The  reason  this  is  allowable  is 
that  the  turbulent  structure  has  an  inherent  self-correcting  character 
which  drives  it  asymptotically  toward  a  limit  which  is  essentially 
independent  of  the  initial  conditions. 

The  persistence  through  space  and  time  of  the  net  overall  influence 
of  the  deterministic  equations  of  motion,  in  addition  to  that  of  the 
continuity  equation,  is  what  imparts  to  the  turbulent  field  its  under¬ 
lying  orderly  structure,  upon  which  are  superimposed  phenomena  of  a  more 
random  character.  One  method  of  expressing  the  orderly  structure  is 

in  terms  of  the  correlation  coefficient  R  defined  in  an  earlier  section. 

•3 

Let  us  here  apply  this  concept  to  a  simple  grid  of  type  A,  with  NJ  grid 
points.  In  the  present  context,  the  appropriate  definition  of  corre¬ 
lation  coefficient  is  in  the  form  of  the  matrix: 


R(p*q>r,T) 


N  N 


N 


7  F  Z  I  |  J  u(i,  j*k,t)*u(i+p,J+q,k+r,t+T)dt 
i«l  j=l  k=l  0 


(12-1) 


& 


(12-2) 


where  p) 

q)  -  0,1, 2, 3, . (N-l) 

r) 

Since  each  of  the  indices  p,q,r  varies  over  N  distinct  values,  the 

expression  (l2-l)  defines  N^  separate  equations.  These  equations  fix 

the  values  of  the  Ir  elements  of  the  correlation  matrix  fl(p,q,r,T). 

Each  element  is  a  function  of  the  time  parameter  f.  This  matrix  of 

functions  provides  a  reasonably  adequate  description  of  the  orderly 

structure  cf  the  turbulent  field.  Ultimately,  the  matrix  functions  can 

be  calculated  from  the  computer  simulation  results. 

Now  consider  the  special  case  T  *  0  which  establishes  all  spacewise 

■5 

correlations  at  any  fixed  time.  The  elements  of  the  matrix  R  are 

■5 

now  fixed  constants.  Suppose  these  NJ  values  are  known  or  assumed,  and 
are  imposed  as  essential  initial  conditions.  This  additional  constraint 
reduces  the  initial  number  of  degrees  of  freedom  from  2N  to  N  .  Pre- 
sumably,  the  remaining  degrees  of  freedom  may  now  be  selected  at 
random  from  some  suitable  probability  distribution  without  appreciably 
affecting  the  structure  of  the  turbulence. 
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13*  PARTIAL  DERIVATIVES  IN  TERMS  OP  SURFACE  INTEGRALS 


In  dealing  with  functions  of  position,  we  usually  need  to  estimate 
not  only  the  values  of  the  function  its._*lf,  but  also  values  of  Its 
various  partial  derivatives.  Now  a  direct  analytical  differentiation 
of  the  approximate  quasi-linear  distribution  as  defined  by  any  of  the 
foregoing  methods,  A,  B  or  C,  is  not  adequate  for  this  purpose.  The 
reason  is  that  such  an  approximate  distribution  is  only  piecewise  con¬ 
tinuous.  Estimates  of  first-order  derivatives  obtained  by  direct 
differentiation  of  this  approximation  would  be  discontinuous;  estimates 
of  higher  order  derivatives  would  be  indeterminate.  This  accords  with 
the  well  known  fact  that  the  derivatives  of  an  approximating  function 
tend  to  be  less  accurate  and  well  behaved  than  the  approximating  func¬ 
tion  itself.  On  the  other  hand,  integrals  of  an  approximating  function 
tend  to  be  more  accurate  and  well  behaved  than  the  approximating  func¬ 
tion  itself.  Fortunately,  it  turns  out  to  be  possible  to  express  the 
derivatives  of  a  function  in  terms  of  certain  surface  integrals,  as 
explained  below.  Estimation  of  derivatives  by  this  method  of  integra¬ 
tion  avoids  the  rapid  deterioration  of  accuracy  associated  with  direct 
analytical  differentiation;  consequently,  the  integration  method  is 
always  used  in  our  work. 

With  the  integral  method  explained  below,  derivatives  are  calcu¬ 
lated  at  all  principal  grid  points;  derivative  values  between  grid 
points  are  then  estimated,  as  usual,  by  quasi-linear  interpolation. 

This  process  may  be  repeated  as  often  as  necessary  to  obtain  estimates 
of  the  higher  order  derivatives. 

In  the  calculus  of  vectors,  the  gradient  of  a  scalar  function  at  a 
point  is  defined  bb  the  resultant  of  a  certain  limiting  process.  First, 
a  small  volume  element  6v  is  defined  whose  surface  ftS  encloses  a  point 
in  question.  See  Fig.  13*1.  Secondly,  an  average  value  of  the  gradient 
within  the  small  enclosed  volume  is  defined.  This  definition  involves 
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(a)  VOLUME  ELEMENT  OF  GENERAL  SHAPE 
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(b)  CUBICAL  VOLUME  ELEMENT 


Fig.  13.1  Volume  Elements  for  Definition  of  Partial  Derivatives 


the  evaluation  of  a  certain  integral  over  the  bounding  surface  of  the 
eie'se-’t.  Finally,  the  volume  of  the  element  is  allowed  to  shrink  to 
zero.  The  resulting  limit  of  the  average  gradient  is  defined  as  the 
local  gradient  at  the  point. 

This  rather  complicated  definition  may  be  written  symbolically  in 
the  form 

ttp  *  lim  t-  f  n/trlS  (13-i) 

tS 

where  dS  is  an  infinitesimal  surface  element  on  the  bounding  surface 
6S,  while  n  is  an  outward  unit  vector  normal  to  area  dS. 

Inasmuch  es  the  gradient  is  expressible  in  terms  of  partial  deri¬ 
vatives  in  the  form 


it  follows  that  the  partial  derivatives  of  cp  are  also  definable  in 
terms  of  surface  integrals,  that  is, 


<i3‘3> 

6s 

<f£>  ■  3  •  -  si"-  0  {  lb  1 3  •  “rfs)  <13J,) 

6S 

<ff>  ■ *  •  ■  ii"-  o  { h  js*  •  W  <i3-5> 

Likewise,  the  divergence  of  a  vector  function  u  may  be  written 
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However,  in  work  with  numerical  approximations  of  the  present  type 
in  which  the  functions  involved  are  specified  only  at  discrete  grid 
points,  and  the  variation  between  these  points  is  merely  estimated  by  a 
kind  of  three-dimensional  quasi-linear  interpolation,  the  above  process 
of  going  to  the  limit  of  zero  volume  cannot  be  accomplished  with  suf¬ 
ficient  accuracy.  As  mentioned  above,  the  previously  described  inter¬ 
polation  schemes  do  not  in  themselves  provide  sufficient  continuity  in 
the  representation  of  the  function  to  define  accurate  local  values  of 
derivatives  at  a  specified  point  by  direct  differentiation. 

Suppose,  however,  that  we  omit  the  above  limiting  process.  Instead 
we  choose  small  but  finite  cubical  elements  as  established  by  the 
finite  spacing  of  the  basic  grid  points.  Now  the  quasi-linear  inter¬ 
polation  scheme  permits  us  to  evaluate  the  various  surface  integrals 
required  in  the  above  definition*  with  excellent  precision.  Hence  we 
may  calculate  very  ac  irate  average  values  of  gradient,  partial  deri¬ 
vative,  velocity  divergence,  etc.,  over  the  volume  element  as  a  whole. 

In  some  instances,  average  values  over  the  element  are  all  that  we 
require.  Thus  in  dealings  with  incompressible  flows,  an  exact  solution 
would  theoretically  require  that  the  velocity  divergence  vanish  at  all 
points  of  the  field,  not  just  at  the  basic  grid  points.  This  ideal  is, 
of  course,  neither  attainable  nor  really  necessary.  For  a  practical 
numerical  solution  it  is  quite  sufficient  to  require  that  the  average 
divergence  vanish  for  each  of  the  small  cubical  cells  into  which  the 
field  is  divided.  Thus  in  Method  A,  for  example,  for  an  overall  cubical 
region  or  block  of  length  L  along  each  edge,  with  a  basic  grid  spacing 
of  dimension  (l/n),  there  are  individual  volume  elements  involved, 
each  a  cube  of  length  (L/N)  on  a  side.  -Hence  the  non-divergence  require¬ 
ment  for  the  field  boils  down  to  requiring  that  the  average  divergence 

3 

vanish  for  each  of  these  N  distinct  volume  elements. 

In  some  cases,  however,  the  concept  of  an  average  value  over  the 
volume  does  not  suffice.  What  is  required  is  a  method  of  establishing 
local  point  values  at  the  basic  grid  points,  so  that  a  complete  spatial 
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distribution  of  the  quantity  may  be  estimated  by  means  of  the  quasi - 
linear  Interpolation  method  previously  mentioned.  For  example,  in 
applying  the  momentum  Integral  equations  to  volume  elements,  ve 
establish  accurate  average  time  rates  of  change  of  the  velocities  over 
these  volumes.  This  is  not  really  sufficient.  What  is  needed,  of 
course,  are  local  values  at  the  basic  grid  points  themselves. 

To  meet  this  requirement,  ve  introduce  one  additional  postulate, 
namely,  that  the  average  value  of  a  quantity  over  a  given  volume 
element  constitutes  the  best  available  estimate  of  its  local  value  at 
the  centroid  of  that  element.  Because  of  this  postulate,  it  is  of 
course  necessary  to  establish  the  pattern  of  control  volumes  and  grid 
points  in  such  a  way  that  the  grid  points  are  in  fact  located  at  the 
centroids  of  the  control  volumes.  When  this  is  done,  th3  formulas  of 
this  section  may  readily  be  applied  to  calculate  the  values  of  diver¬ 
gence,  gradient  and  partial  derivatives  at  all  the  grid  points. 

The  results  obtained  in  this  way  will  depend  on  the  particular 
differencing  and  interpolation  scheme  adopted,  such  as  the  methods  A, 

B  or  C  described  elsewhere  in  this  report.  Detailed  results  for  differ¬ 
encing  method  C  are  presented  in  Appendix  A. 

The  following  results  are  presented  as  typical  examples  of  specific 
differencing  formulas  which  are  obtainable  by  use  of  the  foregoing 
methods.  We  show  the  typical  partial  derivative  (g^)  as  defined  by 
each  of  the  differencing  techniques  A,  B  and  C  explained  earlier.  Other 
partial  derivatives  can  also  be  found  merely  by  systematic  permutation 
of  indices  in  the  following  formulas. 


Method  A 


[<P(i>J>k)  j  *  (|)  {  Ig  [_<p(i+l> Jik)-cp(i-l,J,k)] 

+  37  [»(i+l,J,k-l)4«(i+l,j-lflt)-Kp(i+l,J+l,k)4o(i+i,j,k+i)  (13.7) 

-cp(i-l,  J,k-l)-<o(i-l,  J-l,k)-p(i-l,J+l,k)-cp(i-l,.j,k+l)"| 

+  [o(i+l,J-l,k-l)+to(i+l,J+l,k-l)-Kp(i+l,J-*.,k+l)4tD(i+i,j+i;k+l) 

-$(i-l,j-l,k-i)-cp(i-l,j+l,k-l)-cp(i-l,j-l,k+l)  -co(i-l,  j+l,k+l)  |  } 

Method  B 

[cpX(i,j,k)  |  *  (|)  {  |  [vI(i+l,J,k)-vI(i-l,j,k)  | 

+  rc°II(i /  J >k) ■«DII(i> J-l>k)-KPII( i >  J>k-l)+©XX(i, J-l,k-l) 

-cpII(i-l>J>k)-CDII(i-l,J-l,k)-coIi(i-l,j/k-l)-cpII(i-l;j-l/k-l)j  }  (13.8) 

fj-  [cpTI(i, J,K) J  =  (|)  |VX(i+l,j,k)-cpn(i-l,J,k)j 

+  Jg  [cDI(i+l,J+l,k+l)-K0I(i+l,j>k+l)+OI(i+l/<)+l,k)-KpI(i+l;JA) 

-  cDI(i,J+l,k+l)-crI(i,J,k+l)-nI(i,J+l,k)-toI(i,J/k)  j  ^ 


(13-9) 


Method  C 


§£  [V(i>J>k)J  * 

(lj|){<pII(i,J,K)VI(i>J-1>k)KpII(1,J,K-1)4<pII(1>J-1>k-1) 

■CpII(i-l>  J,k)-cpII(i-l,J-l,k)-cpII(i-l> J,k-l)-cpII(i-l, j-l,k-l)l  (i3_10) 

§£  [cpn(i,J,k)  J  = 

(^){coI(i+1,J+1,k+1)V(i+1,J,k+1)-KnI(1+1,J+1,k)-KDI(i+1,j,k) 

-w^i^+i^+D-co^i^^+DVd^+i^)-®1^^^)}  (13.n) 

It  1b  also  of  Interest  to  consider  a  simple  differencing  formula 
such  as  might  have  been  used  in  the  absence  of  the  finite-element 
principles  expounded  in  this  report.  Thus,  for  a  simple  grid  of  the 
same  type  as  used  in  Method  A,  ve  might  have  simply  written,  by 
inspection. 


5J  [cp(i,J,lO]  -  (|j0 {<p(i+l, J,k) -cp(i-l,  J,k)^  (13-12) 

The  relative  crudeness  of this  approximation  becomes  quite  apparent 
when  this  expression  is  compared  with  the  previous  differencing 
formulas. 
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1U .  CONTINUITY  EQUATION 


A 

Let  U  represent  the  mean-flow  velocity  vector,  and  u  the  turbulent 
perturbation.  We  apply  the  law  of  conservation  of  mass  to  a  small  but 
finite  control  cell  of  volume  tv  and  surface  6S.  Let  n  be  the  outward 
normal  unit  vector  at  surface  element  dS.  For  the  combined  mean  flow 
plus  turbulence,  this  gives 


A  A  p  A  A  A 

V* (U+u)=  I  (U+u) *ndS  =  0  (lU-1) 

t  S 

For  reasons  explained  earlier,  we  have  omitted  the  process  of  going  to 
the  limit  of  zero  volume. 

By  averaging  this  equation  over  time  (or  over  a  statistical 
ensemble  of  macroscopically  similar  cases),  we  find  that  the  mean  flow 
itself  must  satisfy  the  following  mass -conservation  requirement,  namely, 

V*U  =  ~  |  U*ndS  *=  0  (lU-2) 

*>S 

Subtracting  (l^-2)  from  (l^-l)  gives  the  net  mass -conservation 
relation  for  the  turbulence  itself,  that  is, 

V.u  =  D  =  ^  [  u*ndS  =  0  (lU-3) 

The  integral  in  ( lU -3 )  is  evaluated  over  the  surface  of  the  ele¬ 
ment.  But  for  the  cubical  control  volumes  here  considered,  the  surface 

consists  of  six  squares,  each  of  area  6S.  Let  the  individual  surfaces 

L  2 

be  represented  by  the  index  k=l,2,3,^,5>6.  Note  that  area  &S  *  (n-) 

L  q  n 

and  volume  tv  *  (-)  .  Therefore  Eq.  (lU-3)  ntay  be  reduced  to  the  form 
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D  ’  ®  l  {  fe  ]  “•"4S}  ’  0  (lk‘k) 

k-1  6s  k 

The  quantity  endoeed  within  braces  represents  the  average  maBs  flow 
per  unit  area  through  the  kth  bounding  surface* 


15.  MOMENTUM  EQUATION 


We  apply  the  momentum  theorem  of  fluid  flow  to  the  same  small  but 
finite  Eulerian  control  cell  used  in  the  previous  section.  The  follow¬ 
ing  nomenclature  is  used  to  denote  the  various  forces  per  unit  area 
(and  per  unit  mass)  acting  on  the  differential  surface  element  dS, 
namely, 

A 

F  =  mean  viscous  force 
n 

A 

f'  ■  turbulent  viscous  force 
n 

a  °  mean  pressure 

cp  =  turbulent  pressure 

For  the  combined  mean  flow  plus  turbulence,  we  obtain  the  equation 


^  f  (U-Ki)dv  -  -  ^  J  (U+u)(U+u)-ndS 

*  h  i  <VVas  -  k  J  ^is 

fts  ss 


(15-1) 


For  the  particular  case  of  the  simple  uniform  shear  flow  mainly 
considered  in  this  report,  the  following  simplifications  apply  to 
Eq.  (15-1),  namely, 


<K> 

(15-2) 

A 

1*  ndS  =  0 

(15-3) 

' - , 

^i> 

3cl 

in 

n 

O 

(15-4) 
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U  =  iOy 


P  «  * 

I  uu- 


ndS  =  0 


(15-5) 

(15-6) 


Furthermore,  the  turbulent  viscous  forces  fn  may  be  expressed 

A  A  A 

specifically  as  f  ,  f  ,  f  ,  namely,  the  forces  on  faces  normal  to  the 
x  y  z 

x,y,z  axes  respectively.  Theee,  in  tarn  may  be  written  in  the  form 

;  .  i  (  TJZ)  +  3  (  'ja.'S*  i  (  (15-7) 


A  A  T  v  A  T  v  A  ^  T  . 

f  =  1  '  +  J  l  -^)  +  k  f 

z  v  o  y  Vo/  vby 


(15-8) 

(15-9) 


The  viscous  stresses  themselves  can  be  expressed  in  terms  of 
velocity  derivatives.  For  an  incompressible  fluid  of  constant  viscosity 


V  p; 


(ID 

(15-10) 

(  v  (p-  + 

V  p  J  v  p  /  >cy  szy 

(15-13) 

^y.-' 

(15-11) 

(  =  v  ( 
v  p  v  p  y  J^z 

(15-lfc) 

vdZ J 

(15-12) 

(  JSC)  „  f  s  V(*X  + 

V  p  ./  V  p  >  -i*x  5y 

(15-15) 

Of  course,  the  various  partial  derivatives  shown  here  must  be  ex¬ 
pressed  in  terms  of  appropriate  surface  integnls,  as  explained  earlier. 

With  the  various  conditions  listed  from  (15-2)  through  (15-15) 
incorporated  or  implied,  the  basic  momentum  equation  (.15-1)  is  reduc¬ 
ible  to  the  form 
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A  A 


h .!  (??)dv  =  - 0  h  J  y(iu+ui)*nds 


As 


1_ 

6v 


uu»ndS 

fc/ 

6s 


1_ 

Av 


f  ds  -  t~ 
n  Av 


CondS 


AS 


is 


(15-16) 


This  is  the  basic  result  required,  although  it  is  still  in  a  some¬ 
what  unwieldy  form.  It  is  convenient  to  introduce  some  auxiliary  nota¬ 
tion  as  follows.  Let 


h  1  <S)dv  ■  -  “t  (15-17) 

Av 

A 

where  u^  is  the  mean  time  -ate  of  change  of  u  over  the  volume  Av,  and  is 
taken  as  the  best  available  estimate  of  the  local  point  value  at  the 
centroid  of  the  element. 

Let  the  remaining  integrals  be  evaluated  first  over  the  six  indi¬ 
vidual  square  bounding  surfaces,  k=l,2,3,^, 5,6.  Thus,  let 


S,_  = 


T,.  = 


-  0  7  1  y(iu+ui) *ndS^ 

AS  1 

f  1  f  A*  ‘  1 

^  os  •  , 
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-  f  }r  I  cindS1 
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(15-18) 

(15-19) 

(15-20; 

(15-21) 
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Furthermore,  let  ub  introduce  the  abbreviation 


A  A  A  A 

«k  -  +  \  +  *k 


(15-22) 


With  this  notation,  the  basic  momentum  equation  (15-16)  reduces  to 
the  relatively  simple  statement 

6 

\  ■  <!>  I  <VV  <»-23> 

k«l 


Dy  way  of  interpretation,  it  might  be  remarked  that  all  quantities 
(15-18)  through  (15-21)  inclusive  have  the  dimensions  and  character  of 
mean  stress  over  the  kth  face,  that  is,  total  force  on  the  face  divided 
by  area.  Furthermore,  is  the  stress  associated  with  the  shear  rate 
0  of  the  mean  flow.  The  turbulent  apparent  stresses  or  momentum  trans- 

A 

port  terms  are  represented  by  T^.  These  latter  terms  are  quadratic  in 
the  velocities  and  are  the  only  non-linear  effects  present.  Purely 

A  A 

.  It  has  been  shown  that  the  H^, 
although  linear,  involve  the  partial  derivatives  of  the  turbulent 
velocity  perturbations,  rather  than  the  velocities  themselves. 

It  is  convenient  in  the  subsequent  analysis  to  lump  the  terms 

A 

together  under  the  separate  symbol  R^.  The  reason  is  that 
these  particular  components  can  be  calculated  directly  from  the  turbu¬ 
lent-velocity  distribution,  without  any  reference  to  pressure  distribu¬ 
tion.  Determination  of  the  pressure  requires  a  separate  and  lengthy 
calculation,  and  the  resulting  average  stress  associated  with  the  pres¬ 
sures  is  therefore  represented  by  the  separate  symbol  P^. 

Turning  to  Eq.  (15-23)  it  is  seen  that  the  terms  have  the  character 
of  net  force  per  unit  mass,  that  is  to  say,  of  acceleration.  It  turns 
out  to  be  advantageous  to  employ  the  additional  terminology  below, 
namely,  6 

F  -  (£)  £  ^  (15-24) 

k-1 
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<wv 


viscous  stresses  are  represented  by 
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Hence  the  basic  momentum  equation  becomes  finally 

*\  A 

ut  -  F  -  tto  (15  -26) 

A 

The  term  p  represents  that  part  of  the  net  accelerating  force  which  can 
be  calculated  directly  fran  the  velocity  distribution  existing  at  a 
particular  instant  of  time.  The  term  (-Vcp)  is  the  net  pressure  force; 
it  can  be  calculated  only  after  the  pressure  distribution  cp  has  been 
found.  The  unknown  pressure  function  co  is  itself  fixed  by  the  velocity 
distribution,  but  the  relation  is  rather  complex.  However,  it  can  be 
solved  by  the  method  outlined  in  the  next  section. 


16.  TURBULENT  PRESSURE  DISTRIBUTION 


The  momentum  equation  (15-26)  has  so  far  been  usee,  only  to  show 

A 

how  the  local  velocity -time  derivative  u.  can  be  calculated  when  the 

A  T> 

direct  velocity-induced  forces  F  and  the  pressure  forces  (-Vrp)  are 

A 

known.  It  has  been  indicated  that  the  forces  F  can  be  calculated  in  a 
relatively  straightforward  manner  when  the  instantaneous  velocity  dis¬ 
tribution  is  known.  The  pressure  forces,  however,  depend  on  the  detailed 
distribution  of  the  turbulent  pressure  ep,  and  this  quantity  is  so  far 
still  an  unknown.  The  purpose  of  this  section  is  to  show  how  cp  itself 
is  determined  from  the  turbulent  velocities. 

To  show  this,  we  need  only  to  epply  the  divergence  operator  term- 
by-term  to  the  momentum  equation  (15-26).  Thus 


7«U.  *=  7«F  -  7*(^p) 

w 


(16-1) 


The  three  terms  in  this  equation  con  be  expressed  in  expanded  form 
as  follows : 


p  if  * 

V*(Vfcp)  «  T^CP  B  i(^p)*ndS 

6S 
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Consequently,  with  these  definitions,  the  basic  equation  (l6-l) 
may  be  rewritten  in  the  concise  form 

72cd  =  7-F  -  (~)  (16-5) 

This  is  the  partial  differential  equation  which  fixes  the  pressure 

A 

distribution  e.  The  quantity  7»F  on  the  right  side  is  a  known  function 
of  the  turbulent  velocity  perturbations.  The  quantity  (^)  calls  for 
some  comment,  however. 

Strictly  speaking,  for  an  incompressible  flow,  the  divergence  D 
should  vanish  identically,  in  which  case  its  time  derivative  will  also 
vanish.  In  actual  computations,  however,  because  of  round-off  or  other 
errors  of  approximation,  the  calculated  divergence  will  usually  differ 
from  zero  by  some  very  small  quantity  *D.  In  order  to  prevent  such 
small  errors  from  accumulating  and  growing,  it  is  desirable  to  assign 
the  value  of  (|^)  in  such  a  way  as  to  tend  to  liquidate  the  error 
already  present.  Thus,  we  wish  to  eliminate  the  current  divergence 
error  rtD  within  some  small  time  interval  *t'.  But  how  shall  this  time 
interval  be  chosen?  If  the  rate  of  error  correction  is  too  small,  the 
divergence  errors  will  tend  to  accumuiate  and  grow,  perhaps  excessively. 
If  the  rate  of  error  correction  is  too  large,  the  associated  pressure 
perturbations  may  become  excessive;  furthermore,  the  divergence  error, 
instead  of  being  eliminated,  may  be  reversed  in  sign  and  possibly  even 
increased  in  magnitude,  thus  leading  to  computational  instability. 
Therefore,  it  is  advisable  to  normalize  6t  '  with  respect  to  appropriate 
reference  parameters  of  the  turbulence.  Cell  size  ^  and  mean  turbu¬ 
lent  energy  E  represent  the  natural  reference  parameters  for  this  case. 
Hence  we  assign 

<K>  ■  -  (ff>  ■  -  °D  E  ^  6D  (l6-6) 

where  is  a  dimensionless  coefficient  whose  optimum  value  can  be 
established  by  numerical  experimentation  on  the  computer. 
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Consequently,  the  basic  equation  for  the  pressure  distribution  now 
becomes 

T2?  m  7*F  +  CD  |  fiD  (±6-7) 

where  the  right  side  is  a  known  function  of  the  space  coordinates. 

Therefore  pressure  op  is  the  only  unknown  in  Eq.  (16-7)  •  For  accur¬ 
ate  work,  it  is  essential  that  the  last  term  in  this  equation  be  of 
▼wry  small  magnitude. 

This  relation  lends  itself  to  solution  for  the  unknown  pressure  of 
op  by  a  numerical  relaxation  technique.  The  essential  principle  of  the 
method  may  be  summarized  as  follows.  The  calculation  starts  with  an 
arbitrary  initial  estimate  or  approximation  for  the  unknown  function  op. 
The  local  values  of  OP  in  the  vicinity  of  the  first  grid  point  are  then 
adjusted  so  that  Eq.  (16-7)  becomes  exactly  satisfied  for  the  first 
control  cell.  Then  the  values  in  the  vicinity  of  the  second  cell  are 
adjusted  in  a  corresponding  fashion.  This  process  is  continued  syste¬ 
matically  to  cover  the  entire  mesh.  However,  each  time  a  local  cell  is 
adjusted,  this  process  slightly  upsets  any  previously  made  adjustments 
in  adjacent  cells.  Nevertheless,  each  step  of  the  computation  produces 
a  net  reduction  in  the  overall  mean  square  error.  Hence  a  complete 
sweep  of  the  mesh  in  the  above  manner  substantially  improves  the 
original  approximation.  Therefore  the  above  procedure  may  be  repeated 
as  often  as  necessary,  until  the  absolute  residual  error  in  Eq.  (16-7) 
is  everywhere  less  than  same  very  small  pre-assigned  limit. 

Calculation  of  the  turbulent  pressure  distribution  ffl  by  the  above 
iteration  method  is  quite  straightforward  in  principle.  The  boundary 
conditions  on  the  vertical  faces  of  the  block  are  simple  because  of 
the  blockwise  periodicity  of  the  solution.  The  boundary  conditions 
on  the  horizontal  faces  are  somewhat  tricky,  however,  owing  to  the 
relative  sliding  between  successive  rows  of  blocks,  as  explained  in 
another  section  of  this  report.  Owing  to  the  number  of  iterations  re¬ 
quired  to  attain  convergence,  the  pressure  calculation  is,  of  course, 
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quite  lengthy.  It  consumes  by  far  the  greatest  proportion  of  the  total 
computing  time. 


17.  BASIC  SURFACE  INTEGRALS  FOR  A  FINITE  ELEMENT: 
METHOD  C 


This  discussion  applies  to  any  scalar  function  of  position,  say, 
for  definiteness,  the  velocity  component  u.  We  choose  differencing 
method  C  for  this  discussion  because,  as  explained  elsewhere,  this  was 
the  method  used  in  the  program  WRBOCODE,  MARK  V.  We  consider  an 
arbitrary  cubical  element,  of  either  family,  and  of  length  2a  «=  L/N. 

The  origin  of  cool'd inates  is  originally  at  the  volume  centroid,  and 
the  dimensionless  coordinates  ?,  T|,  C  take  on  values  of  +  1  along  the 
bounding  surfaces. 

The  basic  inte rpolation  formula  for  this  case  is 

u  -  A0+A1?+A2nfA3C+A4TK+A5C?+A6?n4-A7?nC+A8(l-§2)(l-T12)(l-C2)  (17-1) 

Consider  the  application  of  this  equation  to  the  bounding  surfaces 
of  the  element.  The  surfaces  in  question  may  be  those  normal  to  either 
the  x,  y  or  z  axes,  respectively;  superscripts  are  used  to  distinguish 
these  three  orientations.  Consequently,  Eq.  (17-1)  may  be  reduced  to 
one  of  the  three  forms  below,  namely, 


aJ  +  A^n  +  A*c  +  a*t£ 

(17-2) 

a£  +  A^C  +  +  A^CS 

(17-3) 

A0  +  Al?  +  +  A3?T1 

(17-JO 

The  four  constants  in  each  equation  depend  only  on  the  surface  in 
question.  They  are  uniquely  determined  by  the  values  of  the  function 
at  the  four  corners  of  the  surface  (this  feature  is  one  of  the  advan¬ 
tages  of  Method  C.) 


Consider  the  specific  example  indicated  in  Fig.  17,1.  We  drop  the 
superscript  notation  for  this  portion  of  the  discussion  since  there  is 
no  ambiguity  in  this  case. 

Upon  substituting  the  coordinates  of  the  points  A,  B,  C,  D  into 
(17-4),  and  solving  the  resulting  four  equations  fcr  the  four  initially 


unknown  constants,  we  obtain 

A1  *  r  [•  UA  +  “B  '  UC  *  uDj  <17‘6> 
A2  *  if  ['  UA  ‘  UB  +  “0  +  “d]  (17'7> 
A3  '  f  [+  UA  -  •  UC  *  “d]  <17'8) 


An  analogous  treatment  may  be  applied  to  the  velocity  components 
v  and  w.  It  i6  con’  enient  to  denote  the  corresponding  constants  by 
B's  and  C's,  respectively. 

Now  a  review  of  the  basic  continuity  and  momentum  equations  reveals 
that  the  surface  integrals  involved  therein  are  of  three  main  types, 
an  example  of  each  type  is  presented  below  for  the  case  of  a  surface 
normal  to  the  z  axis.  The  extension  of  these  results  to  other  velocity 
components  and  other  surfaces  is  easily  made. 

+1  +1  +1  +1 

1  g  [  l  u(ad?)(edri)  =  J  J  j  1  A^A^+A^A^pld^dTi  -  Aq  (17-9) 

-i  .1  -1  -i 

+1  +1  +1  +1 

j  j  ^(arl)(ad?Hadrl>  -  J  J  j  [aq+a^+a^a  ?Tilm?dp  =  |(|)a2  (17-10) 
<2a>  -1  -1  -1  -1 
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+1 

J  uv(ad?)(adTl) 
-1 


+1  +1 

f  1 1  [V*  1?+A2nfA3?Tl][B0+B1?+B2nfB3?Tl]d§dri 


-1  -1 


■  W  5<*lVA2B2>  +  5  *3b3  (17-u) 


Note  that  the  fora  (17-9)  represents  a  simple  linear  average  value 
of  the  function  u  over  the  surface.  All  constants  except  A^  vanish 
from  the  result.  In  turn  A^  is  merely  the  arithmetic  average  of  the 
four  corner  values.  Thi6  result  is  consistent  with  elementary  intuitive 
expectations . 

In  result  (17-10)  all  terms  vanish  except  the  coefficient  of  r|. 

The  method  of  subscripting  here  used  is  such  that  for  a  z  surface  the 
non-vanishing  component  is  A2;  for  an  x  surface,  it  would  be  constant 

Al* 

In  this  connection  y  surfaces  are  exceptional,  for  T|  1b  a  constant 
in  this  case  and  may  be  moved  outside  the  integral,  thereby  reducing 
the  integral  from  the  form  (17-10)  to  the  form  (17-9).  As  a  matter  of 
fact,  it  is  advantageous  in  dealing  with  surface  integrals  to  shift  the 
origin  of  coordinates  irom  the  centroid  of  the  volume  to  the  centroid 
of  the  surface  directly  involved.  In  this  way  we  obtain  'or  any  y  sur¬ 
face  the  result  t|  =  0,  whereupon  the  integral  (17-10)  simply  vanishes. 

Perhaps  the  most  striking  result  is  the  integral  (17-11).  This 
goes  well  beyond  what  could  be  achieved  by  mere  intuition,  or  by  treat¬ 
ment  of  the  surface  as  if  it  were  an  infinitesimal  quantity.  Note  that 
all  cross  product  terms  drop  out  of  the  final  result  but  that  neverthe¬ 
less,  all  eight  original  constants  are  significant  and  are  retained. 

This  undoubtedly  represents  a  significant  contribution  to  the  accuracy 
of  the  met tod,  especially  for  relatively  coarse  meshes.  Recall  that 
the  integral  (17-11)  occurs  in  connection  with  the  turbulent  momentum 
stresses,  and  that  the  accurate  evaluation  of  these  stresses  is  one  of 
the  important  goals  of  this  study. 
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The  results  (17-9) >  (17-10)  aid  (17-11)  may  be  generalized  by 
simple  permutation  of  variables  and  subscripts  so  as  to  produce  all  the 
surface  integrals  required  in  this  analysis,  without  further  recourse 
to  detailed  integration.  In  this  way  the  fundamental  continuity,  momen¬ 
tum  and  pressure  relations  may  be  reduced  to  algebraic  difference 
equations  which  may  then  be  directly  programmed  for  the  computer.  Hie 
essential  features  of  this  algebraic  development,  and  of  the  results, 
are  summarized  in  Appendix  A. 
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18.  SCALE  OF  TURBULENCE 


One  of  the  fundamental  aspects  of  turbulence  is  the  linear  scale 
vhich  characterizes  the  phenomenon.  It  has  already  been  pointed  out 
that  our  method  of  analysis  implies  that,  in  some  sense,  cell  size  must 
be  small,  whereas  block  size  must  be  large  in  relation  to  the  scale  of 
turbulence . 

But  how  shall  scale  of  turbulence  be  measured?  Several  approaches 
have  already  been  indicated.  One  is  based  on  the  concept  of  correla¬ 
tion  distance,  another  on  the  wavelength  spectrum  of  the  turbulence. 

A  third  method  previously  mentioned  is  associated  with  the  idea  of  a 
critical  cell  size  below  which  turbulence  cannot  be  self-sustaining. 

All  of  these  methods,  however,  involve  lengthy  and  intricate  calcula¬ 
tions. 

There  is  clearly  a  requirement  for  a  much  simpler  measure  of  scale, 
which  can  be  routinely  calculated.  To  answer  this  need,  a  concept  based 
on  vorticity  has  been  adopted. 

From  the  point  of  view  of  the  finite -element  method,  the  vorticity 
vector  at  a  grid  point  may  be  defined  by  the  relation 

U)  =  t-  i  nxudS  sec  (l8-l) 

0  v  «J 

fts 

Considering  the  separate  scalar  components  of  vorticity,  we  have 


2 


uu 


2  2  2 

UU  +  uu  +  uu 
x  y  z 


-2 

sec 


Similarly,  for  the  velocity  components  at  any  point 


(18-2) 


V2 


2  2  2 
u  +  v  +  w 


2E 


p  p 

ft  /sec 


(18-3) 


where  E  represents  kinetic  energy  per  unit  mass. 
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Both  these  quantities  may  be  averaged  by  integration  over  the 
entire  volume  of  the  block,  thus  yielding  root-mean-square  values 
y  ^2  and  J2E  respectively. 

Nov  a  definition  of  scale  X  may  be  written  in  the  form 


(10-4) 


X  =  CJ  2 

»  U) 


Note  that  the  radical  has  the  required  physical  dimension  of  length.  A 
dimensionless  normalizing  constant  C  is  introduced  to  permit  adjustment 
of  the  units  of  scale.  The  value  of  C  is  assigned  in  such  a  way  that 
the  magnitude  of  X  acquires  a  clear  physical  significance  for  a  certain 
simple  limiting  case. 

The  case  in  question  may  be  reprssented  by  a  hypothetical  velocity 
distribution  of  the  form  illustrated  (for  v)in  Fig.  10.1.  The  equations 


u  *  P  sin 


sln^ 


(18-5) 


_  .  2«x  .  2nz 

Q  sin  — j—  sin 


(10-6) 


R  sin  sin 


(18-7) 


This  is  seen  to  represent  a  simple  reference  sine  wave  of  wavelength  t. 
The  constant  C  is  now  chosen  in  such  a  way  that  the  scale  parameter  X 
for  thi«  simple  case  turns  out  to  be  exactly  equal  to  the  wavelength  l. 
This  same  value  of  C  is  then  retained  as  a  fixed  constant  for  all  other 
cases  as  well.  Consequently  X  may  be  interpreted  as  a  kind  of  general¬ 
ized  mean  wavelength  of  the  turbulence.  (See  Appendix  C,  Section  4.) 

It  is  now  proper  to  say  concerning  cell  size  normalized  scale 
of  turbulence  X,  and  block  size  L  that  we  should  ideally  like  to  have 

|  «  X  «  L  (10-0) 


00 


However,  if  N  be  sharply  limited  by  computer  memory  capacity,  we  are 
willing  to  accept  the  more  modest  statement 


19*  MEAN  EFFECTIVE  TUREUIENCE  STRESSES 


One  of  the  important  objectives  of  the  present  research  is  to 
improve  our  understanding  of  the  so-called  Reynolds  stresses,  that  is, 
the  mean  effective  values  of  apparent  stress  associated  vith  the  turbu¬ 
lent  momentum  transport.  It  is  well  known  that  turbulent  fluctuations 
can  be  admitted  into  the  Navier-Stokes  equations  of  motion,  and  that  the 
resulting  expressions  can  then  be  averaged  over  time  (or  over  an  ensemble 
of  macros copically  similar  cases).  When  this  process  is  carried  out, 
the  equations  which  are  recovered  are  exactly  similar  to  the  ordinary 
Navier-Stokes  equations  for  laminar  flow,  except  for  the  presence  of 
the  additional  Reynolds  stresses.  In  effect,  this  result  tells  us  how 
the  mean  motion  is  affected  by  the  mean  effective  turbulence  stresses. 
However,  this  is  only  part  of  the  story.  What  is  required  to  complete 
the  solution  is  a  knowledge  of  how  the  effective  stresses,  in  turn, 
are  affected  by  the  mean  motion. 

Sometimes  an  attempt  is  made  to  express  this  latter  relation  by 
means  of  an  analogy  with  the  viscous  stresses.  The  distortional  com¬ 
ponents  of  the  viscous  stress  tensor  are  known  to  be  proportional  to 
the  corresponding  components  of  the  strain  rate  tensor  based  on  the 
mean  motion.  The  constant  of  proportionality  is,  of  course,  the 
ordinary  molecular  viscosity.  It  is  sometimes  assumed,  therefore, 
that  distortional  components  of  the  Reynolds  stress  may  likewise  be 
found  by  multiplication  of  corresponding  components  of  the  mean  strain 
rate  tensor  by  a  suitable  scalar  factor  of  proportionality,  the  so- 
called  eddy  viscosity.  If  this  analogy  were  strictly  valid,  it  would 
represent  a  great  simplification  of  the  problem.  Unfortunately,  it  is 
at  best  only  a  rough  kind  of  approximation.  The  melancholy  fact  seems 
to  be  that  there  does  not  necessarily  exist  any  single  common  scalar 
factor  of  proportionality  between  the  various  Reynolds  stress  compon¬ 
ents  and  the  corresponding  mean  strain  rate  components.  Consequently, 
a  more  fundamental  approach  to  the  problem  is  needed. 
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Since  stress  and  strain  rate  are  tensor  concepts,  it  Is  advantageous 
for  the  present  discussion  to  adopt  Cartesian  tensor  notation  and  con¬ 
ventions.  In  particular,  the  coordinates  x,  y,  z  are  replaced  by  x^, 

Xg>  x^  respectively.  Hence  the  general  component  of  the  Reynolds  stress 
may  be  written  in  the  symbolic  form. 

(iK)“  ■ 


where  i)  ,  p  , 

j)  1,2,3 

Contracting  indices  (and  applying  the  usual  summation  convention),  then 
dividing  by  3  gives 


I  (uiV 


2 

1 


(19-2) 


where  p  represents  the  equivalent  "hydrostatic"  pressure  associated 
with  the  turbulence,  while  E  is  the  corresponding  mean  kinetic  energy 
of  turbulence r 

The  distortional  Reynolds  stresses  are  therefore  expressible  in 
the  form 

(  T1)  -  8U  (  P>  -  VV6ij!e  <19-3> 


where  ■  +  1  if  i  *  J 
«  0  if  i  +  J 


(l9-*0 


We  now  tackle  the  problem  of  identifying  the  specific  parameters 
upon  which  the  values  of  the  stress  tensor  might  depend.  For  this 
purpose  it  is  useful  to  distinguish  four  types  of  turbulent  flows  as 
indicated  in  Table  19*1. 


TABLE  19.1 


lypes  of  Turbulent  Flows 


lype 

I 

II 

III 

IV 

Definition 

Stationary? 

No 

Yes 

Yes 

No 

Isotropic? 

Yes 

No 

No 

No 

Homogeneous? 

Yes 

Yes 

No 

No 

Stress  and  Force 

Characteristics 

of  Mean  Flow 

Distortional 

Stresses? 

None 

Present 

Present 

Present 

Net  Reynolds 

Forces? 

None 

None 

Present 

Present 

Discussion  elsewhere  in  this  report  shows  that  for  iype-I  flows, 
the  distortional  Reynolds  stress  tensor  vanishes  identically.  Or  to 
put  the  matter  more  simply,  there  can  be  no  mean  effective  shear 
stresses  in  an  isotropic  flow  field.  Hence  the  Reynolds  stress  is  of 
the  purely  "hydrostatic"  type  for  this  case.  It  follows  that  type-I 
turbulence  is  too  specialized  in  nature  to  provide  information  regard¬ 
ing  a  general  state  of  stres3. 

Flows  of  type  II,  on  the  other  hand,  are  much  more  interesting  in 
that  non-vanishing  shear  stresses  are  now  definitely  involved.  In 
fact  this  category  represents  the  simplest  possible  case  compatible 
with  the  existence  of  a  general  state  of  stress  at  a  point.  Therefore, 
it  would  seem  to  be  the  most  profitable  case  to  study  initially  for 
the  purpose  of  establishing  the  parameters  which  determine  the  stresses 
at  a  point.  The  simple  shear  flow  which  is  the  subject  of  most  of  this 
report  is  a  particular  example  of  a  Type -II  flow. 
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Of  some  interest*  besides  the  Reynolds  stresses  themselves*  are  the 
net  resultant  forces  per  unit  volume  exerted  by  the  Reynolds  stresses 
upon  a  small  volume  element  enclosing  an  arbitrary  point.  These 
resultant  forces  are  proportional  to  the  gradients  of  the  Reynolds 
stresses.  However*  in  any  homogeneous  flow  field*  the  stress  gradients 
are  everywhere  zero,  of  course,  and  the  net  Reynolds  forces  vanish. 

(in  fact  the  net  viscous  forces  associated  with  the  mean  flow  also 
vanish  in  this  case.)  Consequently*  in  any  homogeneous  flow*  the 
Reynolds  stresses  are  determined  by  the  mean  flow*  but  the  mean  flow 
is  not  influenced  by  the  Reynolds  stresses.  This  fact  represents  an 
important  simplification  of  the  overall  problem. 

Once  the  stress  laws  applicable  to  Type-II  flows  have  been  worked 
out*  the  study  may  be  extended  to  the  more  complex  flows  of  Type  III 
and  IV.  There  appear  to  be  some  grounds  for  the  hope  that  the  laws 
ultimately  established  for  Type -II  flows  will  need  only  minor  elabora¬ 
tion  to  serve  adequately  for  the  more  complex  cases.  In  fact  it  might 
well  be  that  under  certain  circumstances,  the  laws  derived  for  iype-II 
flows  may  be  applied  unchanged  to  the  other  cases;  however,  much  more 
work  needs  to  be  done  before  such  a  conclusion  can  be  asserted  with  any 
degree  of  confidence. 

Another  important  application  of  the  foregoing  stress  laws  arises 
in  the  following  connection.  In  principle*  the  method  of  analyzing 
turbulence  described  in  this  report  depends  for  its  success  on  the 
creation  of  a  mathematical  model  having  a  very  large  number  of  me6h 
points.  Since  not  all  possible  difficulties  can  be  faced  and  conquered 
simultaneously,  it  is  tacitly  assumed  in  much  of  the  discussion  in 
this  report  that  the  number  of  mesh  points  required  for  sa;  .sfactory 
results  lies  within  the  capabilities  of  tic  best  modern  computers,  or 
at  least  within  the  capabilities  of  machines  that  will  become  available 
within  a  few  years.  In  truth,  however*  this  is  probably  an  over- 
optimistic  view.  Hence  there  exists  a  need  for  simplifying  the  method 
to  the  extent  that  it  does  truly  fell  vithin  the  capabilities  of  now- 
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existing  eqaipnent.  Naturally,  any  such  simplification  entails  a  cer¬ 
tain  loss  of  information.  However,  if  the  simplification  is  carried 
out  judiciously,  the  loss  of  information  need  not  be  too  serious. 

The  basis  of  simplification  is  obviously  to  ise  a  coarser  mesh 
than  that  demanded  by  the  "exact1'  theory.  This  means  that  components 
of  the  turbulence  of  wavelength  smaller  than  the  minimum  attainable 
cell  size  cannot  be  explicitly  resolved  in  the  analysis.  However,  if 
the  mean  effective  Reynolds  stresses  associated  with  this  range  of 
wavelengths  can  be  correctly  simulated,  the  resulting  error  in  the 
overall  motion  should  be  quite  negligible. 

As  has  been  discussed  earlier,  It  is  believed  that  the  above 
objective  can  be  accomplished  by  means  of  an  analysis  in  two  stages. 
The  first  deals  with  the  small-scale  turbulence,  the  second  with  the 
large-scale  turbulence.  The  cell  size  of  the  large-scale  motion  is 
taken  as  the  block  size  of  the  small-scale  motion.  In  the  initial 
small-scale  analyses,  the  full  memory  capacity  of  the  computer  is  de¬ 
voted  entirely  to  resolving  the  small -wavelength  turbulence.  From 
these  results,  it  should  be  possible  to  infer  simplified  rules  for 
establishing  the  corresponding  small-scale  mean  effective  Reynolds 
stresses.  Such  stresses  can  subsequently  be  incorporated  into  the 
large-scale  analysis  to  produce  the  required  overall  result. 

For  the  remainder  of  this  section,  we  revert  now  specifically  to 
flows  of  iype  II.  Since  the  Reynolds  stresses  are  uniquely  determined 
whenever  the  mean  motion  is  specified,  we  may  hypothesize  that  these 
stresses  are,  in  general,  some  unknown  functions  of  the  mean  velocity 
components  and  of  their  various  partial  derivatives  of  all  orders. 

Also,  we  observe  that  the  Navier-Stokes  equations  of  motion  are  valid 
with  respect  to  any  inertial  frame.  Furthermore,  any  coordinate  frame 
moving  at  any  constant  velocity  with  respect  to  an  initial  frame  is 
also  an  inertial  frame.  It  follows,  therefore,  that  the  Reynolds 
stresses  cannot  depervi  on  the  velocity  components  themselves  since 
these  are  different  in  different  inertial  frames,  but  must  be  functions 
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of  the  velocity  derivatives  alone.  The  velocity  derivatives  are  in¬ 
variant  with  respect  to  a  shift  of  inertial  frames.  However,  for  any 
homogeneous  flow  field,  all  velocity  derivatives  vanish  except  those 
of  first  order.  We  conclude,  therefore,  that  the  six  Reynolds  stresses 

at  any  point  are  each  functions  of  the  nine  quantities  '  _ i  J  where 

i) 


J) 


1,2,3- 


au 


axi 


/  \ 

The  nine  components  of  the  unsymmetrical  tensor  '  r— “■  ’  may  be 
reorganized  into  six  components  of  the  symmetrical  strain^rate  tensor 
7 j.  and  three  components  of  the  rotation  vector  as  follows.  Let 


7ij  2  L  ax1  3XjJ 


(19-5) 


k  *  1,2,3 


iffi. 


.3xJ 


SXjJ  2  ijk  v  Sx^ 


(19-6) 


where  C  ^  *  +1  if  i,  j,  k  are  in  cyclic  order 

■  -1  if  i,  j,  k  are  in  anti-cyclic  order  (19-7) 

=  0  if  any  two  indices  are  equal 
Fbr  any  given  state  of  the  mean  motion  as  defined  by  the  mean 
velocity  vectors 

Ui  "  Ui  (*1>  *2*  x3^  (19-8) 


the  actual  values  of  the  six  strain  rate  components  7^  and  the  three 
rotation  components  will  depend,  of  course,  on  the  arbitrarily 
chosen  orientation  of  the  axes;  consequently  only  three  of  the  six  com¬ 
ponents  /.  characterize  the  actual  nature  of  the  mean  strain  rate. 

*  J 

This  can  be  readily  seen  if  the  principal  axes  are  chosen  as  coordinate 

axes,  for  in  this  case  the  strain  rate  tensor  reduces  merely  to  the 

•  •  • 

three  principal  strain  rates  7^,  722f  733*  which  now  fuHy  define  the 
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character  of  the  field.  For  the  present  discussion,  however,  it  is  more 
convenient  to  choose  the  axes  in  a  somewhat  different  way;  clearly,  the 
choice  of  a  particular  axis  orientation  entails  no  loss  of  generality  in 
the  phenomenon  itself  -  this  still  represents  a  general  state  of  strain 
rate  at  a  point. 

We  choose  the  axes  x^,  Xg,  x^  in  such  a  way  that  x^  and  Xg  lie  in 
that  principal  plane  which  contains  both  the  algebraically  largest  and 
algebraically  smallest  principal  strain  rates.  However,  axes  x^  and  Xg 
are  chosen  to  be  at  ljh°  with  respect  to  the  principal  strain  rate  axes. 

Of  course,  x^  is  the  third  principal  axis,  the  one  associated  with  an 
intermediate  value  of  strain  rate.  The  situation  is  shown  in  Fig.  19.1* 
Also  shown  are  the  familiar  Mohr's  circle  representations  for  this  strain 
rate  condition,  and  for  the  associated  Reynolds  stresses. 

In  this  particular  frame  of  reference,  since  x_  is  a  principal 
axis,  we  have  7^  =  7^  =  0*  The  non-vanishing  components  are  7^,  T ^ 

7 and  712*  However,  in  this  case  7^  *  7  ^  60  that  a8ain  there  are 
Just  three  independent  degrees  of  freedom  in  the  strain  rates.  Further¬ 
more,  we  choose  the  axes  such  that  7^  i6  positive.  Since  axes  x^  and 
Xg  are  at  4 5°  to  the  principal  axes,  it  follows  that  the  magnitude  of 
7  represents  the  mflXimura  possible  shearing  rate.  Reference  to  Fig. 
19.1,  the  Mohr's  circle  diagram,  should  help  make  this  clear. 

In  discussing  either  the  stress  or  the  strain  rate  tensors,  it  is 
frequently  useful  to  distinguish  between  the  isotropic  and  distortional 
components.  The  distortional  components  of  strain  rate  are  defined  by 


where  7  -  j  7^ 


1  r 

3  ^711 


'a  E  '  -*ij7 


+  722  +  733) 


(19-9) 

(19-10) 


However,  for  the  case  of  an  incompressible  fluid, \hich  i6  the  specific 
case  to  which  this  entire  report  is  restricted,  we  have  zero  divergence, 
that  is 


89 


NROL 102  66 


(a)  Reference  Axes  X)(  X~,  X3  and  Principal  Axes  Xp,,  Xpg.  Xpj 
for  Mean  Strain  Rate  Tensor 


(b)  Mohr's  Circles  for  Mean  Flow  Strain  Rates 


(c)  Mohr's  Circles  for  Reynolds  Stresses 


Fig.  19*1  Reference  Axes  and  Mohr's  Circles  of  Strain  Rate  and  Stress. 

(a)  Reference  Axes  x-p  x3  and  Principal  Axes  xpp  Xpp  Xp^  for 
Mean  Strain  Rate  Tensor. 

(b)  Mohr's  Circles  for  Mean  Flow  Strain  Rates. 

(c)  Mohr's  Circles  for  Reynolds  Stresses. 
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7±i  -  37  -  (7n  +  722  +  733)  ■=  0  (l9-H) 

whereupon  the  distortional  component  becomes  identical  with  the  strain 
rate  itself,  that  is, 

h}  ■  hi 

By  reat^n  of  the  zero  divergence,  which  establishes  another  rela¬ 
tion  of  constraint  among  the  strain  rate  components,  the  number  of  inde¬ 
pendent  degrees  of  freedom  in  the  strain  rate  tensor  is  reduced  from 
three  to  two.  We  may  therefore  at  this  point  rewrite  the  strain  rate 


components  in  the  optional 

alternative 

form 

•  • 

’ll" 

(19-13) 

• 

7  “ 

'23 

0 

(19-16) 

*22*  « 

(19-14) 

731 

0 

(19-17) 

*33  "  2i 

(19-15) 

7 12  * 

ft 

2 

(19-18) 

where  c  and  0  now  represent  the  two  independent  degrees  of  freedom  in  an 
explicit  way.  If  we  set  e  =  0  in  this  result,  we  obtain  the  strain  rate 
tensor  associated  with  simple  two-dimensional  shear  flow,  with  shearing 
rate  0;  in  fact  the  simple  parallel  shear  flow  which  is  the  subject  of 
most  of  this  report  is  a  particular  flow  of  this  type.  The  foregoing 
strain  rate  pattern  with  e  /  0  is  therefore  a  generalization  to  three 
dimensions  of  the  simple  two-dimensional  shear  flow:  the  meaning  of  Cl 
has  been  appropriately  generalized  in  a  corresponding  manner. 

We  may  summarize  the  argument  to  this  point  by  stating  that  for 
homogeneous  flows,  the  Reynolds  stress  components  must  be  functions  of 
the  form 

=-u1uJ  =  f  (c,  Cl,  w1,  (J>2,  o»3)  (19-19) 
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This  relation  may  be  simplified  even  further  by  dimensional  analy¬ 
sis.  Using  reference  parameters  p,  v  and  Cl,  we  set  up  the  dimensionless 
forma  of  the  various  variables  as  follows : 


nr-M 

(19-20) 

•  r 

e#-s 

(19-23) 

(>».£.  1 

(19-21) 

n» .  i 

(19-24) 

Ui 

U*  -  - 

1  ATT 

(19-22) 

lb*  ®i 
i  «  fl" 

(19-25) 

Consequently  we  obtain  at  last 

«$,  wg,  ou*)  (19-26) 

Each  of  the  six  dimensionless  stresses  is  therefore  a  function  of  Just 
four  variables,  as  shown.  This  is  the  general  form  of  the  relation 
which  we  initially  desired  to  at  certain. 

In  the  case  of  a  mean  flow  which  is  purely  two-dimensional,  this 
relation  simplifies  drastically,  for  in  this  case 

e*  -  0,  and  -  iu*  =  0  (19-27) 

Consequently 

T*J  -  -  uj"u*  -  f*  (in*)  (19-28) 

where,  from  considerations  of  symmetry,  two  of  the  six  stresses  vanish, 
namely, 

▼g  -  t|~  -  0  (19-29) 

Hence  for  turbulence  in  a  two-dimensional  homogeneous  mean  flow,  the 
four  dimensionless  non-vanishing  Reynolds  stresses  are  functions  of  only 


a  single  parameter,  iu^,  the  dimensionless  rotation  vector.  This  analy¬ 
sis  makes  evident  the  desirability  of  an  early  investigation  of  the 
two-dimensional  homogeneous  flow  in  which  tv*  varies  over  a  range  of 
values. 

Furthermore,  for  the  special  case  of  the  simple  parallel  shear 
flow  considered  in  this  report 

-  -  |  (19-30) 

Consequently  for  this  particular  case,  the  dimensionless  Reynolds  stres¬ 
ses  reduce  to  four  specific  constants. 

Note  that  the  foregoing  method  is  free  of  any  arbitrary  or  specu¬ 
lative  assumptions  regarding  eddy  viscosity,  mixing  lengths  or  the  like, 
although  it  is,  of  course,  restricted  to  flows  of  Type  II. 

An  interesting  question  arises  in  connection  with  the  principal 
axes  of  stress  anr1  strain  rate.  The  principal  axes  of  the  strain  rate 
tensor  are,  by  definition,  axes  of  zero  strain  rate.  Similarly,  the 
principal  axes  of  the  Reynolds  stress  tensor  are  axes  of  zero  shear 
stress.  According  to  our  usual  rheological  notions  of  stress  and  strain 
we  would  expect  that  along  axes  of  zero  shear  strain  rate  the  shear 
stress  would  also  be  zero.  In  other  words,  the  principal  axes  of  mean 
flow  strain  rate  and  Reynolds  stress  would  be  expected  to  coincide.  But 
do  these  particular  rheological  relations  necessarily  apply  to  the  case 
of  fluid  turbulence?  There  appears  to  be  no  definite  proof  that  they 
do  apply  in  general.  This  is  an  item  on  which  it  should  be  possible  to 
obtain  useful  information  by  the  methods  of  numerical  simulation  con¬ 
sidered  in  this  report. 

Elementary  theory  can  also  shed  some  light  on  this  question.  For- 
this  purpose,  let  indices  i,  j,  k  now  refer  specifically  to  the  princi¬ 
pal  axes  of  strain  rate.  Now  if  any  principal  plane  .  is  also  a  plan 
of  symmetry  of  the  mean  flow,  then  it  follows  from  this  symmetry  that 
t*  =  0  and  t*  =  0.  However,  in  any  such  plane  of  yrnnetry,  it  again  fblkvs 
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that  tu*  *  0  and  urt  =  0.  It  would  appear  that  the  converse  must  also  be 

u  _  _____ 

true.  That  is,  if  'Ju*  *  0  snd  ur*  =  0,  then  also  B  o  and  =  o, 
where  the  three  indices  i,  j,  k  are  all  distinct.  Thus, 

If  ■  0  and  (b£  =  0,  then  -  0  and  t*3  =  0  (19-31) 

if  iU*  ■  0  and  uu*  =  0,  then  T^g  *  0  and  T*^  =  0  (19-32) 

if  uu*  =  0  and  uri*  =  0,  then  =  0  and  t*^  =  0  (19-33) 

It  follows  also  that  if  the  mean  flow  be  irrotational,  the  principal 
axes  of  stress  and  strain  rate  are  then  coincident.  The  implication  is 
that  for  highly  rotational  mean  flows,  the  principal  axes  of  stress  and 
strain  may  cease  to  coincide. 
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20.  C0MRJTATI0NAL  VERSUS  HYDRODYNAMIC  INSTABILITY 


The  discussion  in  earlier  sections  should  make  it  abundantly  clear 
that  considerable  ingenuity  and  care  are  required  to  devise  grid  patterns 
and  calculation  procedures  which  will  yield  maximum  accuracy  and  resolu¬ 
tion  for  a  given  computational  effort  and  cost.  A  special  difficulty 
that  presents  itself  in  various  guises  concerns  the  stability  and  con¬ 
vergence  of  the  calculation  procedure.  This  is  one  of  the  fundamental 
problems  generally  associated  with  numerical  methods.  The  particular 
difficulty  in  the  present  application  is  that  turbulence  itself  is  an 
instability  phenomenon.  There  is  always  the  lurking  danger  of  introduc¬ 
ing  inadvertently  same  suboxc  cause  of  instability  in  the  calculation 
procedure,  and  confounding  the  results  with  the  true  physical  instability 
which  is  turbulence. 

One  common  if  indirect  method  of  coping  with  questions  of  computa¬ 
tional  accuracy  and  stability  is  to  compare  results  obtained  by  numeri¬ 
cal  methods  with  some  known  exact  theoretical  results.  Unfortunately, 
in  the  turbulence  problem,  there  are  no  known  theoretical  results  to 
serve  as  standards  of  comparison. 

An  alternative  is  to  compare  with  experimental  data.  There  is  an 
enormous  amount  of  experimental  data  available  on  various  aspects  of 
turbulence.  Nevertheless  adequate  experimental  data  specifically  appli¬ 
cable  to  the  present  boundary  conditions  might  not  be  available  or 
readily  found.  It  may  eventually  prove  necessary  to  undertake  an  inde¬ 
pendent  experimental  program  for  this  express  purpose. 

A  more  direct  approach  is  to  attempt  to  analyze  theoretically  the 
stability  characteristics  of  the  numerical  method  being  used.  Some 
tentative  guide  lines  are  available  in  this  connection,  although  there 
does  not  seem  to  be  available  any  standard  and  definitive  technique. 

The  present  problem  is  especially  difficult  owing  to  the  great  length 
and  complexity  of  the  basic  difference  equations  involved,  and  owing  to 
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their  non-linearity.  Hence  the  theoretical  approach  to  stability  is  a 
full” fledged  research  problem  in  its  own  right. 

It  is,  however,  a  very  important  problem,  and  efforts  will  be  made 
to  pursue  it  as  circumstances  may  permit.  A  successful  theoretical 
analysis  would  provide  definite  and  rational  criteria  of  stability  and 
should  thereby  improve  the  efficiency  of  the  calculation  procedure  and 
the  reliability  of  the  results  obtained. 

Meanwhile,  the  difficulty  is  being  faced  in  a  pragmatic  way.  Various 
common  sense  safeguards  and  criteria  of  a  heuristic  kind  have  already 
been  provided  in  the  program  HJRBO CODE,  MARK  V.  Thus,  for  example,  the 
degree  of  convergence  demanded  in  the  iteration  procedure  which  fixes 
the  pressure  distribution  is  readily  adjustable.  A  parameter  is  also 
provided  to  limit  the  maximum  displacement  of  any  fluid  particle  in  ary- 
time  step  to  any  desired  small  fraction  of  a  single  cell  length.  Pre¬ 
liminary  numerical  experiments  are  producing  useful  guide  lines  concern¬ 
ing  the  optimum  settings  for  these  parameters,  and  concerning  the  corres¬ 
ponding  level  of  residual  error  or  noise  produced  in  the  calculation 
sequence.  Furthermore,  it  is  well  known  that  the  sovereign  remedy  for 
most  instability  problems  is  a  reduction  in  the  size  of  space  and  time 
increments.  Hence  it  is  believed  that  stability  problems  are  not  likely 
to  defeat  the  goals  of  the  investigation,  although  they  might  very  well 
increase  the  effort  and  cost  required  to  attain  these  goals. 

Incidentally,  it  is  believed  that  the  comparatively  sophisticated 
finite-element  techniques  being  used  in  this  problem  are  far  less  suscep¬ 
tible  to  instability  problems  than  the  looser  methods  based  merely  on 
replacing  differentials  by  finite  differences.  It  must  be  conceded, 
however,  that  no  method  of  differencing,  however  refined,  is  altogether 
exempt  from  the  hazard  of  possible  computational  instability. 
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21.  GENERATION  VERSUS  ANALYSIS  OF  TURBULENCE  DATA 


The  numerical  description  of  a  flow  field  involves  the  use  of  a 
suitable  grid  of  coordinate  points  in  that  field.  For  the  ordinary  case 
of  a  cubical  region  of  length  L  on  a  side,  subdivided  into  cubical  cells 
of  length  L/N  on  a  side  according  to  Method  A,  there  are  N°  cells. 

3 

Since  a  grid  point  is  placed  at  the  centroid  of  each  cell,  there  are  N 

grid  points.  In  the  program  IURB0C0DE,  MARK  V,  a  staggered  scheme  of 

3 

cells  is  employed,  Method  C,  which  involves  2N  independent  grid  points. 
Since  MARK  V  may  use  values  up  to  N  =  8,  the  number  of  grid  points  may 
reach  1024  in  all. 

The  state  of  the  turbulent  flow  at  any  instant  of  time  may  be  ade¬ 
quately  described  by  specification  of  the  values  of  the  three  velocity 
components  and  of  the  pressure  perturbation  at  each  point.  Hence  we 
have,  in  MARK  V,  at  least  8n^  items  of  data  involved  for  each  time  step, 
that  is,  up  to  4096  items  per  step.  (Actually,  certain  additional  items 
of  information  are  also  needed  near  the  top  and  bottom  surfaces,  but 
this  is  immaterial  for  the  present  discussion.) 

Now  the  basic  purpose  of  TURBOCODE,  MARK  V  is  to  follow  the  evolution 
of  the  turbulence  resulting  from  an  arbitrary  initial  perturbation. 

Since  the  number  of  time  steps  involved  in  any  adequate  overall  time 
period  is  large,  it  is  clear  that  we  are  confronted  with  an  enormous 
volume  of  data  -  a  veritable  torrent.  It  becomes  a  major  problem  as  to 
what  to  do  with  this  data,  how  to  store  it,  how  to  process  it,  and  so 
on. 

The  basic  tenet  has  been  adopted  that  the  generation  of  the  basic 
data  and  its  subsequent  processing  for  various  analytical  purposes 
should  not  be  confounded  within  the  same  program.  Accordingly,  TURBOCODE, 
MARK  V  is  intended  solely  to  generate  the  data.  The  processing  involved 
therein  is  minimal  and  is  confined  to  those  few  statistical  features 
which  are  clearly  indispensable  for  providing  a  general  indication  of 


97 


\ 


the  state  of  turbulence.  These  include  (for  each  time  step)  space  mean 
values  of  kinetic  energy,  turbulent  shear  stress  and  scale  of  turbulence. 
Consequently  the  generation  of  the  data  is  accomplished  with  maximum 
efficiency  and  in  minimum  time. 

The  results  thereby  generated  are  recorded  and  stored  on  magnetic 
tape  for  subsequent  analysis  and  processing  by  a  separate  program  or 
series  of  programs,  according  to  the  type  of  information  desired.  How¬ 
ever,  for  conservation  of  tape  and  record it.K  time,  provision  is  made  to 
sample  the  generated  data  and  to  retain  ot\ly  such  selected  portions  as 
are  required  for  the  subsequent  analysis.  This  avoids  the  indiscriminate 
taping  of  huge  amounts  of  data  far  in  excess  of  any  realistic  need. 

In  a  sense,  the  generation  and  sampling  of  data  in  the  present 
study  plays  a  role  comparable  to  that  of  physical  experimentation  in 
conventional  research.  The  subsequent  analysis  and  manipulation  of  the 
data  so  obtained  is  a  more  or  less  distinct  and  separate  step. 
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PARTIAL  SUMMARY  OP  PRINCIPAL  CONCEPTS 


The  mean  flow  field  is  a  steady,  incompressible,  parallel  and  uni¬ 
form  shear  flow  in  which  the  mean  velocity  U  is  of  the  form 

A  A 

U  «  i  fi  y  (22-1) 

A 

where  0  represents  the  constant  shear  rate,  and  i  is  a  unit  vector  in 
the  x  direction.  Physical  boundaries  of  the  flow  are  at  y  *  +  •. 

The  turbulence  is  homogeneous,  that  is,  all  statistical  properties 
of  the  turbulence  are  uniform  over  the  flow  field.  However,  the  field 
is  anisotropic,  that  is,  properties  in  one  direction  are  in  general 
different  from  those  in  another  direction. 

All  physical  quantities  are  non-dimensionalized  in  terms  of  the 
reference  parameters:  density  p,  kinematic  viscosity  v,  and  shear  rate 
fi.  All  possible  uniform  shear  flows  of  the  above  type,  when  non-dimen¬ 
sionalized  in  this  way,  reduce  to  a  single  definite  and  unique  case. 

The  energy  spectrum  of  the  turbulence  is  largely  confined  within  a 
limited  raige  of  wavelengths.  If  we  neglect  some  arbitrary  but  very 
small  fraction  of  the  total  turbulent  energy  at  each  end  of  the  spectrum, 

two  corresponding  limits  of  wavelength  are  thereby  defined,  say  L  . 

mm 

and  Lmwy,  between  which  nearly  all  of  the  turbulent  energy  lies. 

The  flow  field  is  subdivided  into  large  cubical  blocks  of  length  L 

on  a  side.  Each  block  is  in  turn  subdivided  into  many  small  cubical 

control  volumes  or  cells  of  length  (l/n)  on  a  side.  The  integer  N  fixes 

the  fineness  of  the  cubical  mesh.  In  the  ideal  case,  it  is  desirable 

that  L  i  L  and  (l/n)  ^  L  ,  .  This  implies  that  N  i  L  / L  .  . 

max  min  max  min 

Since  the  computer  memory  requirements  are  roughly  proportional  to 
■3 

N  ,  it  will  probably  be  impossible  to  meet  the  above  criterion.  In  that 
case,  the  turbulence  must  be  subdivided  into  small-scale  and  large-scale 
components,  and  analyzed  in  two  successive  stages.  The  method  lies 
beyond  the  scope  of  this  summary,  but  is  detailed  in  the  main  text  of  the 
report. 
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The  large  cubical  blocks  into  which  the  flow  field  is  subdivided 
are  arranged  in  horizontal  rows.  The  horizontal  mid -plane  of  each  row 
has  zero  velocity  with  respect  to  the  mean  flow.  Because  of  the  shear 
rate  fl  of  the  mean  flow,  there  is  a  relative  velocity  of  magnitude  OL 
between  any  two  successive  rows.  At  times  t  ■  0,  1  /fl,  2/0,  3 /0. ... 
all  rows  are  aligned  in  an  unstaggered  cubical  array,  such  that  all 
vertical  faces  of  all  blocks  coincide  and  form  continuous  vertical  planes. 
At  all  other  times  the  blocks  in  successive  rows  are  staggered  with 
respect  to  one  another.  See  Pig.  5* 1* 

At  time  t  m  0,  a  set  of  non-divergent  velocity  perturbations  is 
assigned  at  the  centroids  of  all  the  control  cells.  The  perturbation  is 
arbitrary  except  that  all  possible  wavelengths  in  the  range  of  interest 
are  present.  The  initial  distribution  of  the  perturbation  velocities  is 
taken  as  identical  for  all  blocks,  that  is  all  cells  in  corresponding 
positions  in  all  blocks  are  assigned  identical  perturbations.  Conse¬ 
quently,  the  subsequent  motion  throughout  all  blocks  remains  identical 
thereafter.  This  spatial  periodicity  of  the  solution  makes  the  subse¬ 
quent  boundary  conditions  along  the  block  surfaces  perfectly  definite 
and  unique. 

The  flow  field  passes  through  an  initial  transient  in  which  the 
spectral  distribution  gradually  changes  from  its  arbitrary  initial  form. 
Ultimately,  some  definite  equilibrium  form  and  amplitude  are  attained, 
and  thereafter,  the  turbulence  remains  in  an  essentially  stationary 
state.  The  resulting  velocity  compo..  ;nts  and  pressures  over  the  entire 
block  may  be  sampled  for  any  desired  number  of  successive  time  steps. 

With  methods  suggested  in  the  main  body  of  this  report,  this  data  may 
be  analyzed  to  establish  any  desired  statistical  properties  of  the  turbu¬ 
lence,  such  as  the  mean  kinetic  energy,  Reynolds  stresses,  scale  of 
turbulence,  energy  spectrum,  various  correlation  coefficients,  etc. 

The  initial  state  of  the  system  is  defined  by  the  distributions  of 
the  three  velocity  components  and  the  pressure.  However,  it  may  be 
shown  that  the  equations  of  motion  and  continuity,  along  with  the 
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boundary  conditions,  suffice  to  determine  the  pressure  distribution 
uniquely  whenever  the  velocities  are  fully  specified.  Hence  the  pres¬ 
sure  distribution  is  e.  derived  function  which  cannot  be  specified  inde¬ 
pendently.  Furthermore,  not  all  of  the  velocity  components  are  independ 
ent,  for  it  ma;  be  shown  from  the  condition  of  continuity  and  from  the 
specified  spacewise  periodicity  of  the  turbulence  that  if  any  two  of 
the  velocity  components  are  specified  throughout  the  grid,  the  third 
component  is  thereby  uniquely  determined.  Hence  there  are  just  two 
arbitrary  initial  degrees  of  freedom  for  each  distinct  grid  point  of 
the  block.  Strictly  speaking,  these  arbitrary  degrees  of  freedom  exist 
only  at  the  initial  time  t  =  0,  for  the  entire  subsequent  motion  is 
wholly  determined  thereafter. 

Determination  of  the  pressure  distribution  which  must  co-exist  with 
any  specified  distribution  of  the  velocity  components  involves  a  lengthy 
process  of  computation  by  successive  approximations  until  the  results 
converge  to  within  the  required  degree  of  accuracy.  This  entire  itera¬ 
tion  process  must  be  repeated  every  time  the  velocities  change,  which  is 
to  say,  for  every  time  step.  Hence  the  pressure  computation  takes  by 
far  the  biggest  part  of  the  total  computation  time. 

The  small  but  finite  cells  are  used  as  fixed  control  volumes  in  an 
Eulerian  sense.  The  continuity  and  momentum  equations  are  applied  to 
these  control  volumes.  This  involves  evaluation  of  the  forces  and  of 
the  mass  and  momentum  fluxes  pertaining  to  the  six  square  bounding  sur¬ 
faces  which  enclose  each  cell.  The  evaluate  *  is  accomplished  by  means 
of  complete  and  detailed  surface  integrals.  Similarly,  average  values 
of  momentum  and  energy  enclosed  within  the  cell  are  evaluated  by  means 
of  detailed  volume  integrals.  This  painstaking  integration  technique 
constitutes  the  heart  of  the  so-called  finite-element  method.  It  con¬ 
trasts  strongly  with  most  ordinary  differencing  methods  which  use  much 
cruder  estimates  of  the  required  forces  and  fluxes.  The  difference  in 
accuracy  is  particularly  important  for  relatively  coarse  meshes,  such  s 
those  which  are  necessarily  involved  in  the  present  study. 
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For  evaluation  of  the  above  surface  and  volume  integrals  with  good 
accuracy,  it  is  essential  to  have  a  suitable  interpolation  rule  which 
specifies  how  the  various  quantities  vary  in  the  space  between  the 
principal  grid  points.  For  this  purpose  a  linear  interpolation  rule  is 
employed.  It  is  postulated  that  values  vary  linearly  along  any  straight 
line  parallel  to  any  one  of  the  three  coordinate  axes.  It  turns  out  that 
this  rule  determines  the  distribution  uniquely  throughout  the  entire 
space  when  values  are  specified  at  the  principal  grid  points. 

In  connection  with  the  evaluation  of  various  volume  Integrals,  it 
is  further  postulated  that  the  volume  average  of  any  quantity  within  any 
cell  also  constitutes  the  best  available  estimate  of  the  local  value  of 
that  quantity  at  the  centroid  of  that  cell. 

The  finite-element  technique  is  a  method  of  approximation  which 
employs  a  minimum  number  of  definite  basic  postulates,  and  proceeds 
thereafter  in  a  completely  determinate  and  unambiguous  manner  consistent 
with  these  postulates.  The  necessity  for  any  subsequent  ad  hoc  assump¬ 
tions  is  eliminated.  The  numerical  accuracy  and  stability  thereby 
attained  is  probably  the  maximum  attainable  in  relation  to  the  coarse¬ 
ness  of  the  grid  employed. 

The  basic  idea  of  the  present  calculation  method  is  that  if  the 
velocity  distribution  is  completely  known  at  some  instant  of  time,  the 
new  velocity  distribution  a  short  time  later  can  consequently  be  found. 
First  of  all,  the  pressure  distribution  corresponding  to  the  initial 
velocity  must  be  calculated  by  the  iteration  method  mentioned  earlier. 
Next,  the  equations  of  motion  are  applied  to  the  control  cells  to  estab¬ 
lish  the  time  rates  of  change  of  the  velocity  components  at  the  centroids 
of  the  cells.  Finally,  the  new  velocities  are  computed. 

The  above  cycle  of  calculations  is  repeated  for  every  time  step. 
Hence  the  motion  of  the  system  may  be  followed  indefinitely,  or  as  long 
as  necessary  to  obtain  the  desired  Information. 

A  particular  feature  of  the  above  calculation  procedure  relates  to 
the  maintenance  of  a  continuously  non-divergent  velocity  distribution. 
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Because  of  round-off  and  other  errors,  the  velocity  distribution  at  the 
beginning  of  any  time  step  may  deviate  very  slightly  from  the  required 
condition  of  non-divergence.  This  is  corrected  in  the  following  way. 

A  very  small  correction  term  is  included  in  the  pressure  distribution. 
This  term  influences  the  motion  during  the  subsequent  time  step  very 
slightly  and  in  such  a  way  as  to  tend  to  restore  zero  divergence.  Hence 
there  is  built  into  the  calculation  sequence  a  divergence  correction 
which  continuously  compensates  for  accumulating  round-off  and  other 
errors . 

An  important  question  in  connection  with  any  numerical  method,  in¬ 
cluding  this  one,  pertains  to  the  stability  of  the  calculation  procedure. 
Owing  to  the  great  complexity  of  the  final  difference  equations  in  the 
present  problem,  this  aspect  has  not  yet  received  adequate  study.  How¬ 
ever,  a  number  of  common-sense  precautions  have  been  taken  to  minimize 
the  danger  of  computational  instability.  The  main  precaution  involves 
limiting  the  time  intervals  to  very  small  values  such  that  the  maximum 
displacement  of  any  fluid  particle  is  restricted  to  some  small  fraction 
of  the  cell  size.  It  is  known  that  reducing  the  time  increment  cures 
most  problems  of  computational  instability.  Furthermore,  the  improved 
accuracy  associated  with  the  finite-element  method  itself  provides  an 
enhanced  resistance  to  instability  difficulties. 

The  basic  assumptions  and  concepts  here  summarized  suffice  to 
establish  a  method  for  the  computer  simulation  of  the  detailed  station¬ 
ary  turbulence  in  a  uniform  shear  flow.  The  data  generated  in  this  way 
is  in  some  respects  comparable  to  that  which  might  in  principle  be 
obtained  by  direct  experimental  measurement.  The  computed  data,  how¬ 
ever,  pertains  to  a  very  fundamental  case  but  one  which  would  be  rather 
difficult  to  set  up  experimentally.  Furthermore,  the  computed  results 
are  incomparably  more  comprehensive  than  any  that  could  reasonably  be 
obtained  by  experiment.  The  detailed  data  so  generated  represents 
fundamental  information  which  may  subsequently  be  analyzed  from  various 
points  of  view  to  establish  corresponding  overall  statistical  and 
phenomenological  characteristics  of  the  turbulence. 
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23.  TYPICAL  OCMFUTER  RESULT 


Table  23. 1  and  Pig.  23. 1  show  sane  results  of  a  typical  computer 
run.  For  Internal  consistency  this  data  Is  Included  merely  as  an 
example;  no  detailed  conclusions  are  offered  at  this  time.  More  compre¬ 
hensive  data  and  analysis  are  planned  for  the  near  future.  Nevertheless, 
the  particular  result  shown  supports  the  Idea,  advanced  on  theoretical 
grounds  In  an  earlier  section  of  this  report,  that  for  block  size  L 
below  sane  critical  value,  Lcr,  the  initial  turbulence  dies  away  and 
laminar  flow  Is  restored.  The  data  also  shows  that  the  normalized 
divergence  errors  (column  DBIG)  and  the  inherent  noise  level  In  the 
computer  solution  are  extremely  small. 

The  data  listed  in  Table  23. 1  represents  only  the  gross  overall 
features  of  the  turbulence.  (Most  of  the  columns  are  self-explanatory: 
ENERGY  «  energy  of  turbulence  per  unit  mass;  UVBAR  «  Reynolds  shear 
stress;  SCALE  ■  scale  of  turbulence  as  discussed  elsewhere  in  this  report; 
NITER  ■  number  of  pressure  Iterations  for  convergence.)  Not  shown,  but 
available  from  this  same  computer  inn,  is  an  enormous  mass  of  detail 
showing  the  individual  velocity  components  and  pressures  at  every  point 
in  the  grid  for  every  time  step. 
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Table  23. 1  Typical  Computer  Output 
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PROBLEM  number  32 
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Fig.  23.1  Dissipation  of  Turbulent  Energy  for  Dub-Critical  Block  Size. 


106 


24.  POSTSCRIPT 


The  purpose  of  this  postscript  is  to  record  two  important  concepts 
whose  value  and  significance  were  not  fully  appreciated  at  the  time  the 
main  body  of  this  report  was  written. 

The  first  idea  pertains  to  a  simplified  analysis  of  turbulence  by 

treatment  on  a  purely  two-dimensional  basis.  One  could  set  w  «=  0,  ~ 

oz 

(any  function)  *=  0,  and  treat  u,  v  and  o  as  functions  of  x,  y  and  t  only. 
Actually  this  approach  was  considered  earlier,  but  was  initially  rejected 
on  the  ground  that  true  turbulence  is  inherently  three  dimensional  and 
should  be  treated  as  such.  Nevertheless,  it  appears  that  an  analysis  of 
two-dimensional  "pseudo  turbulence"  would  be  of  considerable  value  and 
interest  as  an  interim  step  toward  the  final  goal.  It  would  undoubtedly 
resemble  true  three-dimensional  turbulence  in  many  significant  respects. 
As  a  two-dimensional  problem  in  the  space  coordinates,  it  falls  well 
within  the  range  of  existing  computer  capabilities.  It  would  enable  us 
to  increase  the  important  mesh  resolution  parameter  N  without  involving 
prohibitively  large  computer  storage  requirements.  It  would  also  facili¬ 
tate  all  subsequent  steps  such  as  spectral  analysis  and  the  like.  While 
obviously  not  a  complete  substitute  for  the  full  three-dimensional 
treatment,  it  would  certainly  be  a  valuable  complement  to  the  latter. 

It  would  be  particularly  useful  and  most  relevant  in  the  earlier  stages 
of  investigation  as  a  tool  for  establishing  basic  concepts  and  methods. 
The  idea  of  using  simplifications  of  this  type  is  nothing  new;  somewhat 
similar  simplifications  have  been  made  by  Bnmons  (Ref.  2),  Kraichman 
(Ref.  3)  and  others. 

The  second  concept  is  a  'further  generalization  of  the  idea  of 
equivalent  reference  frames,  as  outlined  in  Section  5  of  the  report. 

In  Section  5  the  basic  unit  of  fluid  to  which  the  concept  of  equivalence 
has  been  applied  has  been  the  cubical  block  of  length  L  on  a  side.  But 
it  is  both  possible  and  advantageous  to  go  further.  We  can  apply  the 


107 


\ 


concept  of  equivalence  to  the  individual  cubical  cell,  of  length  ^  on  a 

N 

side.  In  this  case,  the  condition  of  equivalence  can  be  reduced  to  the 
statement  that  the  centroid  of  each  cell  shall  have  zero  velocity  with 
respect  to  the  mean  flow.  This  amounts  to  adoption  of  a  kind  of 
Lagrangian  System,  in  the  6ense  explained  in  Appendix  D.  The  scheme 
has  obvious  merit  in  that  all  control  volumes  are  truly  equivalent.  The 
equations  of  motion  are  identical  in  form  for  all  control  cells;  there 
are  no  terms  needed  to  reflect  differences  in  mean  flow  velocity  at  the 
centroids.  Consequently,  the  complexities  and  peculiarities  of  the  slid¬ 
ing  boundary  condition,  which  is  associated  with  the  Eulerian  system  of 
rigid  reference  blocks,  simply  vanish  from  the  analysis. 

If  the  control  cells  are  rigid  cubes  of  length  (^)  on  a  side,  and 
if  the  centroids  of  these  cells  move  with  the  main  flow,  the  effect  is 
as  shown  in  Fig.  26.1.  At  the  initiel  (dimensionless)  time  t*  «=  0,  and 

at  subsequent  integral  values  of  t*  ■  1,  2,  3,  . ,  the  cells  form  a 

simple  cubical  array  as  shown.  At  other  times  the  individual  layers 
become  staggered  as  shown.  The  basic  unit  of  analysis  is  now  not  a 
rigid  block,  but  a  stack  of  slabs  measuring  L  by  L  horizontally,  by 
(jp  vertically. 

The  differencing  formulas  to  be  used  in  this  case  will  contain  the 
time  variable  t,  since  they  must  account  for  the  continuously  changing 
geometry  of  the  system.  However,  the  basic  method  of  deriving  these 
formulas  in  terms  of  surface  integrals  still  applies  in  the  usual  way. 
There  seem  to  be  no  basic  difficulties  involved.  Nevertheless,  some 
care  is  needed  in  defining  the  interpolation  spaces  which  now  also  vary 
with  time. 

Consider  the  set  of  perallelopipeds  formed  by  connecting  principal 
grid  points  by  straight  lines.  In  horizontal  planes  these  grid  lines 
form  squares  of  size  in  the  usual  way.  In  vertical  planes  z  -  constant 

however,  owing  to  the  continuous  shearing  motion,  the  "vertical"  lines 
of  the  grid  lie  at  an  angle  a  with  respect  to  the  y  direction,  as  shown. 
This  angle  changes  continuously  with  time.  We  impose  the  restriction. 
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Fig.  26.1  Differencing  Method  for  Grid  Points  Which  Move  with  Mean  Flow. 
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however,  that  at  all  times  the  angle  a  must  lie  between  the  limits 

-  1/2  s  tan  a  s  +  l/2  (26-1) 

These  limits  are  reached  at  dimensionless  times  t*  =  l/2,  3/2,  5/2,  . 

Each  time  the  upper  limit  tan  a  =  +  l/2  is  reached,  the  lines  are  all 
instantaneously  flipped  back  to  tan  Ct  *  -  l/2  in  order  to  satisfy  (26-1). 
The  angle  then  smoothly  increases  again  to  the  upper  limit,  whereupon 
the  above  instantaneous  flip-flop  action  is  repeated.  The  corresponding 
grid  lines  define  a  system  of  time -varying  interpolation  spaces.  The 
moving  lines  form  the  edges  of  these  spaces.  It  is  quite  natural  to 
adopt  the  rule  that  within  each  interpolation  space,  the  variation  of 
all  quantities  shall  be  linear  along  any  line  parallel  to  any  grid  line. 

The  use  of  interpolation  spaces  whose  shape  varies  with  time  in  a 
linear/ cyclic  pattern  suggests  that  the  control  volumes  might  be  allowed 
to  change  in  shape  in  an  exactly  analogous  manner.  This  is  indeed  en¬ 
tirely  feasible,  and  probably  advantageous  from  a  computational  stand¬ 
point.  At  any  rate,  the  choice  between  rigid  or  time -dependent  control 
cell  ::hape8  can  be  made  on  the  basis  of  convenience,  since  either  choice 
is  correct.  However,  time-dependent  shape  of  control  cells  seems 
slightly  more  consistent  with  the  time-dependent  shape  of  the  interpola¬ 
tion  spaces.  If  this  4dea  is  adopted  consistently,  we  will  have  at 

times  t*  *  0,  1,  2,  .  all  control  cells  and  interpolation  spaces 

instantaneously  cubical.  Similarly,  at  times  t*  «  l/2,  3/2,  5/2,  ...., 
all  control  cells  and  interpolation  spaces  flip-flop  instantaneously 
from  tan  a  ■>  +  l/2  to  tan  cc  *  -  l/2.  However,  the  principal  grid  points 
themselves  continue  to  move  on  with  a  majestic  Lagrangian  dignity, 
utterly  unperturbed  by  the  periodic  flip-flop  of  the  reference  grid 
lines. 

The  gist  of  this  postscript  is  that  there  is  an  early  requirement 
for  a  simplified  two-dimensional  treatment  of  turbulence,  and  that 
future  work  should  employ  coordinates  which  move  with  the  mean  flow. 
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APPENDIX  A 


DEVELOPMENT  OF  BASIC  DIFFERENCE  EQUATIONS 
(METHOD  C) 


A.l  LOCATION  OF  GRID  POINTS 


In  each  of  the  horizontal  directions  x  and  z,  and  for  both  families 
of  points,  the  number  of  grid  rows  needed  is  N.  In  the  vertical  or  y 
direction,  however,  two  extra  rows  are  needed  at  the  bottom  and  two  at 
the  top,  to  meet  the  special  boundary  conditions  at  these  locations. 
Hence  the  number  of  stations  needed  in  the  vertical  direction  is  (N+M. 
Consequently  the  indices  i,  j,  k  in  the  two  families  vary  over  the 
following  ranges 

i1  =  i11  =  1,  2,  3, . N 

J1  «  J11  '  1,  2,  3, .  (N+M  (Al-1) 

k1  =  k11  =  1,  2,  3, . N 


Spacing  between  any  two  successive  rows  in  the  same  family  is  (L/n). 
Any  point  (ijJjk)1*  of  family  II  is  staggered  by  the  amount  (L/2N)  in 
each  of  the  positive  x,y,z  directions  with  respect  to  the  corresponding 
point  (i,J,k)*  of  family  I.  In  the  vertical  direction,  the  lowest  row 
is  at  J1  =  1,  the  highest  at  j1*  =  (N+*0;  the  reference  plane  y  =  0  is 
chosen  to  lie  exactly  midway  between  the  lowest  and  highest  grid  rows, 
for  the  sake  of  symmetry  in  the  S  forces.  In  the  x  and  z  directions, 
the  origin  of  coordinates  is  immaterial,  and  is  arbitrarily  placed  at 


=  0. 


The  following  grid  point  locations  satisfy  the  above  requirements: 
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J 

I 


N 


N 


N 


<*-!> 

II 

Xi  = 

h 

II 

l  .  ») 

i  2 

(k-  |) 

h 

(Al-2) 


The  constants  C*  *  ^  and  ^  in  the  above  results  are  found 

from  the  conditions 


y11  -  y1  =  -  -  (c11 
yj  yJ  N 


Cl>  "  +  k 


(Al-3) 


I  II 

'l  ' 


*  *  ^  *  +  I  <5  -  c*  ’  c“>  '  0 


■II, 


Locations,  other  than  those  at  the  principal  grid  points  themselves, 
are  denoted  by  dimensionless  coordinates  ?,  rj,  C.  Thus  within  a  volume 
cell  of  family  I,  with  centroid  at  (i,  J,  k)1,  or  anywhere  along  the 
bounding  surfaces  of  thi6  cell,  which  are  surfaces  passing  through  points 
of  family  II,  we  have 

-1  s  5  ■  h  (x  ■  xi)  4  +1 

-l  5  *1  -  (y  -  yj)  5  +i 

-1  s  «  ‘  5L  4  41 


Similarly,  within  a  volume  cell  of  family  II,  with  centroid  at 
(i,  J,  k)1  ,  or  anywhere  along  the  bounding  surfaces  of  this  cell,  which 
are  surfaces  passing  through  points  of  family  I,  we  have 
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1  s  ?  -  §l  (xj1  -  x)  s  +1  (Al-5) 

1  s  r\  -  ~  (y*1  -  y)  <  +1 

1  s  c  ■  Ie  <*?  ■  *> s  -11 
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A.  2  SURFACE  CONSTANTS 


Surfaces  passing  through  points  of  family  I,  and  perpendicular  to 
the  x,  y  or  z  axes,  are  denoted  by  superscripts  lx,  Iy  or  Iz,  respect¬ 
ively.  For  surfaces  passing  through  points  of  family  II,  we  use  super¬ 
scripts  I  lx,  Ily  or  Hz  in  like  manner. 

The  three  -/elocity  components  u,  v,  w  are  assumed  to  vary  linearly 
along  any  line  parallel  to  the  x,  y  or  z  axes.  Hence  along  a  typical 
surface  -  let  us  say  a  surface  lx  -  we  have 

u  -  aJx  ♦  Af  n  ♦  a*  c  +  a-  nr 

v  =  bJx  +  n  +  B 2X  c  +  b*x  nr  (A2-1) 

w  «=  cjx  +  c*x  n  +  c*x  C  +  c*x  nc 

where  the  A's,  B's  and  C's  are  constants  whose  values  are  to  be  deter¬ 
mined  from  ohe  known  values  of  the  u's,  v's  and  w's  at  the  four  grid 
points  which  constitute  the  corners  of  the  surface.  (This  approach, 
Method  C,  ignores  the  additional  data  available  from  the  two  other  nearby 
grid  points,  lying  off  the  surface  on  either  side  of  its  centroid.) 

Equations  similar  to  (A2-1)  may  also  be  written  for  th^;  other  five 
surfaces  denoted  by  superscripts  Iy,  Iz,  IIx,  Ily,  IIz. 

Taking  for  illustration  only  the  expression  for  u  in  (a2-1),  we 
have  at  the  four  corners 

u  =  u1  (i,J,k)  n  =  0 

u  *  u*  (i,J+l,k)  n  -  +1 

u  =  u*  (i,J,k+l)  n  ■  0 

u  =  u^  (i,J+l,k+l)  n  = 


C  -  o 

C  -  0 

C  »  +1 
C  =  +1 


(A2-2) 
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On  the  other  hand,  for  the  corresponding  surface  of  family  II,  we 


have  at  the  corners 

u  *=  u11  (i,J,k) 
u  =  u11  (i,J-l,k) 

r\  ■  +1 

T1  =  0 

C  =  +1 

C  *  +1 

(A2-3) 

u  =  u11  (i,J,k-l)  13  ®  +1  C  =  0 

u  *=  u11  (i, J-l,k-l)  r|  =  0  C  0  0 

Now  from  (A2-1)  and  (A2-2)  the  solution  for  the  A's  becomes 

AjX(i,J,k)  -  £  [uI(i,j,k)+uI(i,j+l,k)+uI(i,J,k+l)+uI(i,J+l,k+l)  i 

A*X(i,J,k)  *  j_-uX(  "  )+uI( 

"  )-uI( 

’’  )+ux(  "  y 

1  (A2-UA) 

A^U^k)  -  l  [-uI(  w  )-uX( 

”  )™x( 

”  )-™r(  ”  >] 

A^(i,J,k)  -  J  [+uI(  "  )-uX( 

*  )-uX( 

"  )+uX(  ”  )] 

1 

Similarly  for  the  other  family 


■?*(i,Jik)  -  J  [+uI(i,J-l,k-l)+uI(i,t),k-l)+uI(i,<)-l,k)+uII(i,J,k)J 

^(i^ik)  «  J  [V(  ” 

)+uJ( 

"  )-uX(  "  )+uX  (  "  )J 

r 

n  (A2-UB) 

'2IX(i>JA)  -  l  L^(  " 

>-uX( 

"  )+uI(  ”  )+uI  (  " 

3IX(i,J,k)  -  J  [+uJ(  " 

)-uX( 

”  )-uX(  ”  WX  (  ’’  )] 

The  above  results  may  be  extended  to  the  other  surfaces  by  cyclic 
permutations  as  symbolized  by  the  diagrams 
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i 


X 


A 


y 


? 


In  other  words,  starting  with  any  valid  equation,  and  changing  each  vari¬ 
able  or  index  thereof  by  a  single  step  forward  in  the  above  diagrams, 
produces  another  valid  equation. 

On  the  other  hand,  the  subscripts  as  in  A q,  a^,  A^,  A^  do  not 
change  during  the  above  permutations. 

Furthermore,  the  results  can  also  be  extended  to  the  other  velocity 
components  by  the  permutations  symbolized  below,  namely 


u 


A 


(A2-6) 


In  sum  we  can  say  that  Eqs.  (A2-1+A)  and  A2-1+B),  along  with  the 
permutation  rule6  (A2-5)  and  (A2-6),  fix  all  the  A's,  B's,  and  C's  on 
the  x,  y  and  z  surfaces  of  both  families.  The  final  expressions  obtained 
are  in  terms  of  the  u's,  v's,  or  w's;  that  is,  in  terms  of  the  velocities 
at  the  corner  grid  pointu. 

In  the  sections  which  follow,  various  quantities  needed  for  the 
continuity  and  momentum  equations  are  derived  and  expressed  in  terms  of 
the  above  surface  constants  -  the  A's,  B's  and  C's.  If  desired,  these 
constants  can  then  be  eliminated  from  the  resulting  expressions  by 
substitution  of  the  equations  of  this  section.  Hence  the  results  re¬ 
quired  can  be  expressed  directly  in  terms  of  the  velocity  components 
themselves.  Such  expanded  formulas  were  actually  used  in  IURBOCODE, 

MARK  V  in  order  to  avoid  the  large  computer  storage  which  would  be  re¬ 
quired  for  storing  the  many  surface  constants,  (it  can  be  shown  that  t 
number  of  storage  places  required  for  this  purpose  is  approximately 
72N^(N+4).)  The  expanded  formulas  are  very  long,  however,  and  it  is 
considered  impractical  to  give  them  in  this  appendix. 
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A. 3  DIVERGENCE 


\ 


In  terms  of  the  surface  constants,  the  expressions  for 
in  the  two  families  become 

DX(i,J,k)  -  |  {  [^(i^k)  -  A^U-l^k)] 

+  [s^i^k)  -  ^(i^-l^)] 

+  [c^Iz(i,J,k)  -  cJIz(i,J,k-l)]  } 

Dn(i,j,k)  -  l  {  [Aj*(i+l,J,k)  -  Aj*(i,J,k)] 

+  [Bq^ (i, J+l,k)  -  ^(i,j,k)] 

+  [cJS(i,J,k+l)  -  cJZ(i,J,k)]1- 


divergence 


(A3-1) 


(A3-2) 
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A. 4  mean  stresses  frcm  imposed  shear  rate 


The  shear  stress  components  associated  with  the  interaction  of 
the  turbulent  eddies  with  the  mean  shear  rate  Q*  *  1,  and  defined  by 
Eq.  (15-18),  are  as  follows.  For  family  _ 

sJX(i,J,k)  -  -  (|)  {(J  -  l  -  §)2AqX  (i,J,k)  +  |  (i,J,k)} 

S^U^k)  »  -  (|)  {(J  -  l  -  1)^(1, J,k)  +  \  B^(i,J,k)}  (AU-1) 

sf(i,J,k)  *=  -  (£)  {(J  -  l  -  1)^(1, J,k)  +  l  C^d^k)} 

S ^(i,J,k)  «  -  (§)  (  J  -  f  -  fjB^d,  j,k) 

S^i^^k)  =  0  (A4-2) 

S^(i,J,k)  «=  0 

S^Z(i,J,k)  -  -  (§)  {(J  -  £  -  |)cJz(i,J,k)  +  |  c£Z(i,J,k)} 

SyZ(i,J,k)  »  0  (A4-3) 

S^Z(i,J,k)  -  0 
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Similarly,  for  family  II 

S^i^k)  *  -  (|)  {(J  -Tf  -  f)  SAj^Ci^k)  Ap(l,.J,k)} 
Sy^l^A)  “  •  (^{(J  -J  •  |)  |  B^l^k)} 

s^d,^)  -  -  (§)  {(j  - 1)  1  cj^d^,*)} 

S^d,J,k)  -  -  (|)  (J  -  ?  -  |) 

sf^did**)  -  0 
¥ 

SIIy(l,J,k)  “  0 

Z 

s^Izd,j,k)  -  -  (|)  {(j  - 1)  c£Izd,j,k)  + 1  cfzd,j,k)} 

SIIz(i, J,k)  -  0 

y 

S^Iz(i,J,k)  -  0 


{Pih-k) 
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A. 5  MOMENTUM  FLUX 


The  components  of  the  quadratic  momentum  flux  stresses  defined  by 
Eq.  (15-19)  are  given  below 

r  2 

T^X(i,j,k)  -  -  [_A(JX(i,J,k)J+  |  [A^(i,J,k)] 

+  5  [^(i^k)]  +  i  [A3X(i'«)>k)]  ^ 

^(ifjfk)  -  -  {B^(i,j,k)  A^(i,J,k)  +  J  B^fi,  J,k)A^ X(i  ,  J  ,k) 

+  |  BgX(i,  J,k)  A*X(i,j,k)  +  |  B^(i,J,k)A*X(i,  j,k)} 

T^(i,J,k)  -  -  {c^d^kjA^d^k)  +  |  (^(ijJ^A^U^k) 

+  J  cJX(ilJ,k)AjX(i,J,k)  +  |  C^X(i,J,k)A^C(i,J,k)}  (A5-1) 


These  results  suffice  to  fix  the  pattern.  The  remaining  equations 
may  be  found  by  application  of  the  permutation  rules  (A2-5)  and  (A2-6). 
Exactly  similar  expressions  apply  to  the  surfaces  of  family  II;  it  is 
necessary  only  to  change  superscript  I  to  II.  There  are  nine  separate 
equations  in  each  family,  or  eighteen  equations  in  all  for  the  T  com¬ 
ponents. 


A. 6  LOCAL  VISCOUS  STRESSES 


The  local  viscous  stresses  at  a  point  follow  the  pattern  below. 

For  family  I,  the  local  normal  and  shear  stresses  are,  respectively, 
of  the  form 

T^(i,J,k)  -  ®  [^(i^k)  -  AjIX(i-l,J,k)]  -  |  DI(i,J,k) 

tJ2(i,j,k)  - 1  {  [cj^d,^)  -  <£*(1,4-1*)] 

(A6-1) 

+  [l£Iz(i„J,k)  -  BjIz(i„J,k-l)]  } 

For  family  II,  the  corresponding  forms  are 

T^(i,J,k)  -  [A^i+l,  j,k)  -  A*X(i„J,k)]  -  |  DIX(i,  J,k) 

Ty2(i,J,k)  »  (|){  [c^d^+l, k)  -  0^(1, J,k)] 

+  [Bjz(i,J,k+l)  -  BqZ(1, J,k)]  }  (A6  ^ 

The  other  four  stresses  of  each  family  are  Jound  by  the  usual  per¬ 
mutations  of  (a6-1)  and  (a6-2)  . 

The  divergences  D*  and  which  occur  in  these  equations  are  of 
the  order  of  the  round-off  error  only  and  may  be  neglected  if  desired. 
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A. 7  MEAN  VISCOUS  STRESSES 


The  mean  effective  viscous  stresses  over  the  surface  as  defined  by 
Eq.  (15-20)  are  computed  from  the  local  viscous  stresses  according  to 
the  following  format 
For  family  I 

H^(i, J,k)  «  J  +  T^(i,J+l,k) 

+  T^(i,J,k+D  +  T^(i,J+l,k+l)] 

I^X(i,J,k)  -  l  [^(i^k)  +  T^(i,J+l,k) 

+  ^(i^k+l)  +  T^(i,j+l,k+l)] 

H^(i,J,k)  -  £  [r^(i,J,k)  +  T^x(i,J+l,k) 

+  T^(i,J,k+l)  +  T^(i,J+l,k+l)]  (A7-1) 

The  components  on  the  other  two  surfaces  are  now  found  by  the 
usual  pe mutation. 

For  family  II  only  the  indices  differ.  The  first  equation  suffices 
to  establish  the  sequence 

+  -  t”(l,J,k)]  (A7-2) 

etc. 
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The  local  stresses  may  be  eliminated  from  (A7-1)  and  (A7-2)  by 
use  of  (A6-1)  and  (a6>2).  The  process  Is  lengthy  but  not  difficult* 

A  single  typical  result  Is  the  equation 

-  S  {AJ7jC(i,Jik)  ♦  aJ^U^IA)  +  AjIx(l,J,k+l) 

+  A^CijJ+ljk+l)  -  *£^(1-1,4,10  -  AjtSc(l-l,J*l,k) 

-  Aj^Cl-l^A+l)  -  i^d-lii+lA+l)} 

■  5  +  DI(i,j+l,k)  +  DI(l,J,k+l) 

♦  DI(i,J+l,k+l)}  (A7-3) 

The  extension  of  these  results  to  all  eighteen  components  follows 
In  the  usual  way. 
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A. 8  RESUI/EANT  PRIMART  STRESSES 


The  sum  of  the  S,  T  and  H  stresses  as  given  by  Eq.  (15-22)  Is 
termed  the  resultant  primary  stress  R.  These  R  stresses  are  of  the 
form 

-  S^U^k)  +  T^(l,^k)  +  1^(1, J,k)  (A8-1) 

Eq.  (A8-1)  represents  the  x  component  on  the  lx  face.  Components  In 
the  y  and  z  directions  are  Indicated  by  corresponding  subscripts,  and 
the  other  surfaces  are  Indicated  by  corresponding  superscripts.  There 
are  eighteen  R  stresses  In  all.  Note,  however,  that  eight  of  the 
eighteen  S  components  In  Eq.  (a8-1)  are  zeros. 
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A. 9  NET  PWMAHX  FORCES 


According  to  Eq.  (15*24)  the  net  primary  forces  are  found  by  sum¬ 
mation  of  the  primary  stresses  over  the  six  bounding  surfaces  of  the 
volume  element.  The  results  are  of  the  following  farm. 

For  volume  elements  with  centroids  of  family  I,  bounded  by  sur¬ 
faces  of  family  II  - 

-  (|)  { 

+  [R^(l,J,k)  -  R^(l,J-l,h)] 

+  [RjIZ(i,J,k)  -  R£IZ(i, j,k-l)]  }  (A9-1) 

Similarly,  for  volume  elements  with  centroids  of  family  II, 
bounded  by  surfaces  of  family  I 

-  (|)  {  [r^U+I,^)  -  R^X(l„J,k)] 

+  [l^diJ+lik)  -  R^(l,J,k)] 

+  [r^Z(1,J^+D  -  R^Z(l,J,k)J  }  (A9-2) 

Of  course  the  corresponding  y  and  z  components  are  found  by  chang¬ 
ing  all  subscripts  accordingly  In  Eqs.  (A9-l)  and  (A9-2).  There  are 
three  force  components  In  each  family,  or  six  In  all. 


126 


A. 10  MEAN  PRIMARY1  FORCES 


Note  from  Eq.  (16-3)  that  in  order  to  determine  the  divergence  of 
the  net  primary  forces,  it  is  necessary  to  evaluate  quantities  of  the 
form 


G  -  J  P*n  dS  (A10-1) 

AS 

We  call  these  quantities  the  mean  primary  forces.  Ry  the  usual  inter¬ 
polation  rules  of  Method  C,  the  mean  value  over  the  surface  6S  is 
simply  taken  as  the  mean  of  the  four  corner  values. 

Thus  for  surfaces  of  family  I,  ve  obtain  the  form 

G?  ‘  J  +  +  +  ^(i,J+l,k+l)]  (A10-2) 

and  similarly  for  the  y  and  z  components. 

Thus  for  surfaces  of  family  II 

0-  -  £  [^X(i,J,k)  +  F“(i,M,k)  +  F^d^k-l) 

*  if  d,J-i,k-i)]  (,1a.,) 


and  similarly  for  the  y  and  z  components. 

Xx  Xy  Xz 

The  total  number  of  components  is  six,  namely,  0  ,  0  ,  0  , 

TTv  TTv  Hz  X  y  Z 

G  ,  0  ,  G  .  Note  that  the  subscript  and  superscript  letters 

X  jT  Z 

are  always  in  agreement. 
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A.  11  PRESSURE  COMPATIBILITY  EQUATION 


It  is  convenient  for  this  section  to  adopt  a  simple  letter  symbol 
for  the  divergence  quantity  V*F.  We  arbitrarily  choose  the  letter  Q 
for  this  purpose. 


Prom  Eq.  (16-3) 

ve  obtain 

QI(i,J,k)  - 

(?) 

*k)  -  0^(i-l,J,k)l 

X  j 

+ 

[(^(i^k)  - 

a^d,j-i,k)] 

+ 

[ofZ(i,J,k)  - 

o”Z(l,J,k-l)]  } 

(All-1) 

QIX(i,J,k)  - 

(|)  {  [o*(i+l,J,k>  -  of (i,J,k)] 

4 

•  [o^(i  ,J+l,k) 

-  o^(i,j,k)] 

4 

•  [o^(i,j,k+l) 

-  G*Z(i,J,k)]  } 

(All-2) 

Substituting  these  results  into  Eq.  (16-7)  then  fixes  the  pressure 
compatibility  equation,  namely 

?V(i,J,k)  -  QI(i»J>k)  +  Cp  |  J~H  DI(i,J,k)  (All-3) 

^cpn(i,J,k)  -  Qn(i,J,k)  +  Cn  rfli  Dn(i,J,k) 


(AH-*) 


These  results  fix  the  required  values  of  T^cp  throughout  most  of 
the  grid.  Near  the  top  and  bottom  of  the  block,  however,  special  con¬ 
siderations  apply  because  of  the  sliding  boundary  conditions.  In  the 
first  place  the  pressure -compatibility  equations  (All -3)  and  (All-4) 
are  not  even  needed  for  the  four  rows  J  ■  1,  J  ■  2,  J  ■  N+3,  J  ■  N+4. 
The  reason  is  that  these  are  extra  rows  outside  the  block.  The  velo¬ 
cities  in  these  rows  are  not  computed  by  the  usual  method,  but  are 
found  indirectly  by  taking  advantage  of  the  sliding  blockvise  periodi¬ 
city  of  the  phenomena.  Hence  values  of  V2^)  are  needed  only  over  the 
range  3  <  J  <  (N+2).  It  turns  out  that  Eq.  (All -3)  is  valid  over  this 
range  except  for  the  single  vtlue  J1  ■  3.  Similarly,  Eq.  (All -4)  is 
valid  over  the  required  range,  except  for  the  single  value  J1*  •*  N+2. 
Modified  equations  are  required  for  these  two  particular  rows.  The 
necessary  special  equations  are  derived  in  Appendix  B. 

Tracing  back  the  Q  terms  in  Eqs.  (All -3)  and  (All -4)  shows  that 
they  are  definite  and  known  functions  of  the  turbulent-velocity  dis¬ 
tributions.  Hence  they  represent  a  constraint  on  such  as  to  make 
the  resulting  distribution  of  the  pressure  function  cp  compatible  with 
the  already  existing  velocity  perturbations.  Hie  last  term  on  the 
right  of  (All -3)  and  (All -4)  represents  a  correction  for  error  diverg¬ 
ence.  The  positive  dimensionless  constant  should  be  large  enough 
to  provide  good  divergence  correction,  but  not  so  high  as  to  cause 
computational  instability;  its  optimum  value  is  best  found  by  numerical 
experimentation  cn  the  computer  Itself. 


A.  12  SOLUTION  FOR  PRESSURE  BY  ITERATION 


By  application  of  the  differencing  formula  (13-10)  we  can  express 
3qp/3x  »  cp^  at  any  grid  point  in  terms  of  the  values  of  cp  at  adjacent 
grid  points. 

2  2 

A  repetition  of  this  procedure  gives  3  cp/dx  *  3cp  /3x  *  cp  at 

X  XX 

any  point  in  terms  of  adjacent  grid  point  values  of  cp.  By  cyclic  per¬ 


mutation  of  subscripts,  corresponding  results  can  be  found  for 
2  2  2  2/ 

3  cp/3y  *  3cp  /dy  «  cp  and  for  3  cp/3z  ■  3cp  /3z  ■=  cp  .  Adding  these 

r  yy  z  ^  zz 

three  results  finally  gives  an  expression  for  v  cp  =  (cp^  +  cp^  +  cpzz) 
at  any  grid  point  in  terms  of  the  values  of  cp  at  nearby  grid  points. 


The  procedure,  while  lengthy,  is  not  difficult. 

It  is  convenient  to  multiply  through  by  the  factor  ^  (l/n)2  in 
order  to  simplify  the  expression.  In  this  way  we  obtain  the  result 


2  p 

I  (|)  *V(i,J,k)  -  flpI(i, J,h)  +  \  jyu-l^k)  +  Cpx(i,j-l,k) 

+  «pI(i,J,k-l)  +  cpI(i+l,J,k)  +  tpX(i,J+l,k) 

+  cpI(i, J,k+l)"j  -  |  [cpX(i-l,J-l,k-l)  +  cpI(i-l,J-l,k+l) 

+  cpI(i-l,J+l,k-l)  +  cpI(i-l,J+l,k+l)  +  cpI(i+l,J-l,k-l) 

+  cpI(i+l,j-l,k+l,  +  <pX(i+l,J+l,k-l)  +  cp1  ( i+1 ,  J  +1 ,  k+l)  J 

-  J2  +  cpI(i-l,j,k-l)  +  cpI(i-l,  J,k+l) 

*  cpI(i-l,J+l,k)  +  cpX(i,  J-l,k-l)  +  cpI(i, J-l,k+l) 

+  cpI(i,J+l,k-l)  +  cpX(i,J+l,k+l)  +  cpI(i+l,J-l,k) 

+  cpI(i+l,j,k-l)  +  cpI(i+l,J,k+l)  +  cpI(i+l,J+l,k)]  (A12_1) 
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Exactly  analogous  results  apply  to  the  other  family;  it  is  neces¬ 
sary  only  to  change  all  superscripts  from  I  to  II.  Hence  a  single 
explanation  will  serve  to  cover  both  cases. 

Suppose  we  have  available  a  set  of  trial  vulues  of  the  function  co 
at  all  grid  points.  If  this  is  in  fact  the  correct  function,  then  it 
will  satisfy  Eq.  (A12-1)  at  all  grid  points;  the  left  side  of  this 
equation  is  known,  of  course,  from  the  calculations  of  the  previous 
section.  However,  if  the  trial  solution  is  not  correct,  the  terms  or. 
the  right  side  of  (A12-1)  will  not  in  general  equal  the  known  values 
required  on  the  left. 

We  must  therefore  devise  a  method  for  adjusting  the  cp's  on  the 
right  so  that  the  required  value  on  the  left  can  always  be  obtained 
at  any  specified  point  (i,J,k).  This  procedure  is  then  applied  in 
succession  to  each  point  in  the  field,  for  both  families  of  points. 

The  entire  field  is  swept  in  this  way  again  and  again  until,  at  length, 
the  values  of  the  required  corrections  become  everywhere  less  than 
some  very  small pre-as signed  error.  This  then  fixes  the  required  pres¬ 
sure  distribution. 

An  essential  constraint  upon  the  solution  is  that  the  average 
pressure  for  each  family  of  points  is  zero.  This  follows  from  the 
very  definition  of  the  turbulent  pressure  perturbation  as  the  net 
deviation  from  the  prevailing  mean  pressure.  Consequently  either  the 
procedure  for  adjusting  the  cp's  should  leave  the  mean  value  unchanged, 
or  it  must  be  followed  eventually  by  another  correction  which  restores 
the  zero  mean.  This  restoration  is  easily  accomplished  by  addition  of 
the  appropriate  constant  to  the  values  of  c o  at  all  points  in  the  field. 
An  additive  constant  leaves  all  differences  unaffected. 

The  simplest  method  of  adjustment  is  to  correct  the  value  of 
cpI(i,J,k),  which  is  the  first  term  on  the  right  side  of  Eq.  (A12-1), 
so  as  to  bring  this  equation  into  balance.  This  is  the  method  that 
has  been  adopted  in  TURBOCODE,  MARK  V.  This  procedure  pulls  the  solu¬ 
tion  away  from  zero  mean,  but  this  is  easily  corrected  later  in  the 
manner  indicated  above. 


131 


it  is  not  difficult  to  modify  the  adjustment  procedure  so  that 
it  does  not  disturb  the  mean  value  in  the  first  place.  For  this 
purpose  it  is  necessary  only  to  adjust  not  only  o*(i,J,k),  but  also 
the  next  group  of  six  poincs  as  veil.  These  six  points  are  all 
adjusted  by  equal  amounts.  This  additional  degree  of  freedom  may  be 
used  to  keep  the  mean  value  unchanged  while  Eq.  (A12-1)  is  brought 
into  balance. 

More  elaborate  relaxation  procedures  are  also  possible.  Note 
that  the  27  terms  on  the  right  side  consist  of  the  term  cp*(i,J,k) 
itself,  plus  three  other  groups  of  six,  eight  and  twelve  members, 
respectively.  If  a  separate  correction  quantity  is  assigned  to  each 
of  these  four  categories,  this  provides  four  degrees  of  freedom.  One 
of  these  is  used  to  satisfy  Eq.  (A12-1),  another  to  maintain  zero  mean. 
The  remaining  two  may  be  utilized  to  minimize  the  square  of  the  error. 
This  leads  to  perfectly  definite  relaxation  formulas.  It  is  somewhat 
dubious,  however,  whether  these  more  elaborate  procedures  have  any  real 
computational  advantages.  Hence  they  have  not  been  utilized  up  to  now. 
However,  the  matter  may  be  worthy  of  further  study. 

Each  time  the  solution  is  advanced  by  one  time  increment,  the 
above  relaxation  procedure  must  be  repeated.  The  previously  existing 
pressure  distribution  is  taken  as  the  first  approximation  for  the  new 
pressure.  If  time  increments  are  small,  the  pressure  distribution  does 
not  change  much  during  any  one  time  step,  and  the  pressure  solution 
converges  in  comparatively  few  iterations. 

Attention  is  invited  to  the  fact  that  the  pressures  cp1  and  cp11 
appear  to  be  virtually  independent  of  each  other.  This  arises  from  the 
manner  of  differencing  in  method  C.  Such  apparent  separation  cannot 
occur  in  differencing  methods  B  or  A.  This  circumstance  might  also 
make  the  pressure  Iteration  procedure  in  method  C  somewhat  more  suscep¬ 
tible  to  instability  difficulties  than  would  be  the  case  with  either  of 
the  other  methods. 
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A.  13  NET  PRESSURE  FORCES 


Once  the  pressure  distribution  cp  has  been  found,  the  net  pressure 
forces  are  easily  determined.  The  relevant  difference  equations  for 
the  two  families  follow  the  patterns  below,  namely 

c£(i,J,k)  -  +  ^cpII(i,J,k)  +  cpn(i,M,k)  +  cpn(i,j,k-l) 

+  cp11  (i,J-l,k-l)  -  cpIX(i-l>  J,k)  -  ®^I(i-l, J-l,k) 

-  oI]:(i-l,J,k-l)  -  pn(i-l,J-l,k-l)}  (A13-1) 

cp^d.^k)  «  ^  ^(i+l^+l^+l)  +  CDI(i+l,J,k+l)  +  coI(i+l,J+l,k) 

+  cpI(i+l,j,k)  -  cpI(i,lj+l,k+l)  -  CDI(i,J,k+l)  -  ©I(i,J+l,k) 

-  ^(i,.}^)}  (A13-2) 

The  other  two  components  in  each  family  are  found  by  the  usual 
cyclic  permutations  of  these  equations. 
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A.14  velocity  time  derivatives 


The  local  velocity  time  derivatives  are  now  of  the  form 

u*(i,J,k)  -  ^(i„J,k)  -  <j£(i„),k) 

Cyclic  permutation  gives  the  corresponding  results  for  v*  and 
The  expressions  for  family  II  are  exactly  similar. 


13^ 


(A14-1) 


A. 15  NEW  VELOCITIES 


The  nev  velocities  at  the  end  of  a  time  increment  T  are  of  the 

form 

uJ(i,J,k)  -  uT(i,J,k)  +  TuJ(i>t),k)  (A15-1) 

The  extension  to  the  other  components  and  to  the  other  family  is  the 
same  as  in  the  previous  section. 

This  step  completes  the  calculation  of  the  new  velocity  distribu¬ 
tions  at  the  end  of  the  time  interval.  The  whole  process  is  then 
repeated  for  the  next  time  interval. 
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APPENDIX  B 


DIVERGENCE  CORRECTION  AT  SLIDING  BOUNDARIES 


B.l  SLIDING  INTERPOLATION 


It  has  already  been  pointed  out  in  Appendix  A  that  in  the  y  direc¬ 
tion  calculation  of  velocities,  forces,  pressures,  and  so  on  proceeds 

in  the  normal  manner  only  over  the  range  j  *  3>  4,  5>  .  N+2.  Values 

outside  this  range  are  obtained  by  an  interpolation  process  which  takes 
into  account  the  sliding  boundary  conditions.  Thus  values  at  j  =  1  and 
J  ■  2  are  obtained  by  this  sliding  interpolation  method  from  the  known 
values  at  J  =  N+l  and  j  *  N+2,  respectively.  Similarly,  the  values  at 
J  «  N+3  and  J  ■  N+4  are  obtained  by  sliding  interpolation  from  the 
known  values  at  J  «  3  and  j  ■  4,  respectively.  It  will  be  shown  that 
this  circumstance  affects  the  divergences  and  the  required  divergence 
corrections  in  rows  J*  ■  3  and  J**  «  N+2. 

The  interpolation  procedure  is  based  on  the  fact  that  the  entire 
solution  is  blockwise  periodic.  This  means,  for  example,  that  at  any 
time  t 


u^(i,2,k)  ■  u1  (i+Nt,N+2,k) 
uIX(i,N+3,k)  «  uIT(i-Nt,3,k) 


(Bl-1) 


The  complication  arises,  however,  that  the  quantities  (i+Nt)  and  (i-Nt) 
are  usually  not  integers,  and  therefore  usually  represent  positions 
which  lie  somewhere  between  grid  points.  Nevertheless,  values  at  such 
intermediate  locations  can  be  estimated  by  linear  interpolation 


136 


between  the  known  values  at  the  two  adjacent  grid  points  on  either  side 
of  the  desired  location. 


Any  real  numbers  can  always  be  divided  into  two  parts,  one  of 
which  is  an  integer,  the  other  a  fraction  of  magnitude  smaller  than 
unity.  It  is  convenient  to  introduce  the  notation 

I(s)  =  integral  part  of  s  (3 

F(b)  =  fractional  part  of  s 

from  which  it  is  obvious  that  for  any  s 

s  ■=  I(s)  +  F(s)  (Bl-3) 

It  follows  from  this  that  the  required  intermediate  point  falls 
between  two  adjac^ut  grid  points  according  to  the  relations 


I  (i  +  Nt)  <  (i  +  Nt)  <  I  (i  +  Nt)  +  1 


or 


I  ( i  -  Nt)  S  (i  -  Nt)  <  I  (i  -  Nt)  +  1 


(Bl-4) 


The  various  integers  in  (Bl-4)  may  or  may  not  fall  within  the 
range  from  1  thru  N,  inclusive.  However,  because  of  the  periodicity 
of  the  solution  we  may  shift  each  integer  by  any  exact  integral  multiple 
of  N  such  that  the  shifted  value  does  lie  within  the  above  range.  To 
express  this  idea  we  introduce  a  shift  operator  S  (  )  such  that 


1  *  S(s)  -  s  +  KN  <  N 


(Bl-5) 


where  K  is  an  integer  whose  value  is  uniquely  fixed  by  the  above 
inequalities. 

With  this  notation  we  may  now  define  the  following  quantities; 


\ 


t'  -  F(t) 

iL  -  S  [l(i+Nt')3 


t *  *  truncated  time 


ijj,  -  S  [l(i+Nt ')  +  l] 
iy  -  S  Cl(i-Nt  ')D  | 

iyp  -  S  [l(i-Nt ')  +  1]  ' 


Shifted  indices  (integers) 

for  interpolation  at 

upper  and  lower  (Bl-6) 

boundaries 


»U  ■  V  ' 

»L  "  “up  ■  1  ”  ?(»t ') 


Weighting  factors 


Consequently,  the  relations  implied  by  (Bl-1)  can  now  be  written 
explicitly  in  the  form 

uX(i,2,k)  -  wLuI(iL^N+2,k)  +  wLpuI(iLp,N+2,k) 

un(i,N+3  |k)  -  WyJ^iy^^)  +  wupun(iup,3,k)  (Bl-7) 

Note  that  these  interpolation  formulas  are  based  on  the  values  of 
the  i’s  and  w's  as  evaluated  at  the  beginning  of  the  time  interval.  The 
subsequent  discussion  will  show  the  usefulness  also  of  interpolation 
formulas  based  on  i's  and  w's  evaluated  at  the  end  of  the  time  step. 

Using  subscript  N  to  denote  new  values  at  the  end  of  the  time 
step  t,  we  write 

tjj  ■  t  +  t  New  time 
t ■  F(t„)  Truncated  new  time 

ijif  S  CKi+NtjJ)] 

S  [l(i+NtjJ)+l] 

■  s 

*UHt  -  s  *  13 


New  shifted  indices 


(Bl-8) 
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New  weighting 
factors 


W  B  V 

UN  LPN 

V  e  V 

LN  UFN 


F(Nt  'N) 
l-F(Nt  'N) 


Of  course,  for  the  functions  being  interpolated,  like  velocities, 
pressures,  forces,  and  so  on,  only  the  initial  values  are  known  at  the 
particular  stage  of  calculation  here  involved.  Hence, for  these  quanti¬ 
ties,  only  the  initial  values  are  involved,  in  any  case.  To  distinguish 
the  following  interpolated  values  from  those  of  (Bl-7)  we  add  the  sub¬ 
script  s.  Then 

u^(l,2,k)  =  w][WuI(iLN,N+2,k)  +  wLpNuI(iLpN,N+2,k) 

Ug1  (i,N+3,k)  =  y\,3,k)  +  wUPNuII(iUPN#3^)  (B1_9) 


In  certain  calculations  it  will  be  useful  to  evaluate  the  small 
differences  between  the  values  given  by  the  two  different  interpolation 
formulas  (Bl-7  and  Bl-9).  Hence  we  introduce  the  further  definitions.. 


AuI(i,2,k)  -  ug(i,2,k)  -  uX(i,2,k) 
Aun(i,N+3,k)  =  ugI(i,N+3,k)  -  uIX(i,N+3,k) 


(Bl-10) 


A  similar  notation  is  applied  to  interpolated  forces  and  pressure 
gradients. 


139 


■KVMt 


B.2  DIVERGENCE  CORRECTIONS 


The  divergence  at  the  lower  boundary  o  »  3  niay  be  written  out  in 
full  as  follows: 

DZ(l,3,k)  - 

+  uII(i>3>k"l)  -  uIZ(i-l,3,k)  -  uXX(i-l,3,k-l) 

+  vXX(  "  )  +  vXX(  "  )  +  vXX(  "  )  +  vXX(  "  ) 

+  wn(  "  )  -  wXX(  "  )  +  wXT(  "  )  -  wXX(  ”  )} 

+  (^;)  {uXX(i,2,k)  +  uXX(i,2,k-l)  -  uXX(i-l,2,k)  -  uXX(i-l,2,k-l) 

-  vn(  "  )  -  vXX(  "  )  -  vXX(  "  )  -  vXX(  "  ) 

+  wXX(  "  )  -  wXX(  "  )  +  wXX(  "  )  -  wXX(  "  )}  (B2-1) 

To  facilitate  subsequent  operations,  we  adopt  the  following  abbrevi¬ 
ated  notation  to  represent  the  above  result,  namely, 

DX(i,3,k)  -  (Jj;)  {uXX(i,3,k)  +  — }  +  (jL)  uXX{(i,2,k)  +  — }  (B2-2) 

This  equation  represents  the  situation  at  sane  arbitrary  time  t, 
regarded  as  the  initial  time  of  the  current  time  step.  Using  subscript 
N  to  denote  new  values  at  the  end  of  the  time  step,  we  may  write 

nj(i*3,k)  -  (Jj)  {uj“(l,3,k)  +  — }  +  (Jj)  {ujx(i,2,k)  +  — }  (B2-3) 


We  subtract  these  two  equations  and  impose  the  condition  that  the 
new  divergence  shall  be  zero.  Then  upon  rearrangement  - 

Dj(i,3,k)  -  0  -  DZ(l,3,k) 

+  {  [u5I(1»3#k)  -  uH(i,3,k)]  + . -} 

(B2-4) 

+  {  [UNI^1,2,k^  "  uII(i>2>k)]  +  . } 

The  velocity  changes  on  the  right  side  of  (B2-4)  are  computed  from 
the  equation  of  motion.  The  situation  is  regular  at  J  ■  3.»  therefore 

[«JI(i,3,k)  -  uIX(i,3A)]  *  t  [r[I(i,3>k)  -  coXI(i,3,k)]  (B2-5) 

However,  at  J  =  2,  the  interpolation  rule  applies.  In  particular, 
the  new  velocities  here  are  of  the  form 

ujr(i,2,k)  - 

-wm  {un(im,N+2,k)  +  T[^I(iLN,N+2,k)-lpXI(iLN,N+2,k)l  } 

*WLFN  (U  ^iLPN,N+2,k^  +  T  (iLFN,N+2'k) 

-  ^I(iLHf'N+2'k)]  }  0*2-6) 

Upon  regrouping  of  the  terms  on  the  right  side  of  this  equation, 
and  employment  of  the  definitions  of  (Bl-9)>  this  result  may  be  re¬ 
written  in  the  concise  form 

uJX(i,2,k)  «  uXI(i,2,k)  +  T  jV^(i,2,k)  -  ©”(i,2,k)  ] 


lUl 


(B2-7) 


\ 


Next  we  subtract  the  original  velocity  to  find  the  net  velocity  change. 
Using  the  notation  of  (Bl-10) ,  we  obtain 

LuJX(l,2,h)  -  un(i,2,k)j 


•  AuJ^i^k)  +  T  [lg(i,2,k)  -  q£Xd,2,k)  -  A«fg(i,2,k)]  (B2-8) 


It  is  instructive  to  compare  this  result  with  (B2-5)*  This  show., 
the  difference  in  the  velocity  change  between  a  regular  interior  point 
J»3,  and  a  point  J»2  exterior  to  the  block. 

Upon  substitution  of  the  expressions  of  type  (B2-5)  and  (B2-9) 
into  (B2-4)  end  rearrangement,  the  basic  divergence  equation  may  be 
reduced  to  the  following  form.  The  various  groups  of  terras  are  identi¬ 
fied  by  number,  to  facilitate  the  subsequent  discussion. 

(fe  { K1*1'3'*) +  — ]  +L®?(l'2'k)  +  ""’I) 

♦  H^l)  ♦ 

-  *Vd!i!k)  -  (Jj;)  {  [^(i^lk)  +  — -] 

+  [^(lAA)  +  — -  ]  }  +  if  ^d^ik) 

+  ~  (j^)  (auII(1AA)  +  ~~  }  (B2-9) 

This  result  establishes  both  the  required  value  of  ^cp,  term  (3), 
and  the  corresponding  pressure  solution.  The  required  value  of  ^cp 
is  fixed  by  the  terms  on  the  right,  which  depend  only  on  the  existing 
velocity  distribution.  Terms  (4)  and  (5)  together  represent  the  main 
effect.  This  is  exactly  of  the  same  form  as  found  in  the  interior 
regions  of  the  block,  except  that  the  term  (5)  involves  the  modified 
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interpolation  formulas  as  signified  by  the  s  subscript.  Term  (6)  is 
the  usual  divergence-correction  term.  The  new  term  is  (7)  which  repre¬ 
sents  an  additional  correction  associated  purely  with  the  sliding  - 
boundary  condition. 

On  the  left  side,  the  term  (l)  is  found  to  represent  the  usual 
solution  for  which  may  be  expressed  directly  in  terms  of  the  pres¬ 
sures  themselves.  Term  (2)  is  a  small  correction  associated  with  the 
sliding  boundary  condition.  This  correction  must  be  included  in  the 
iteration  solution  for  the  pressures,  in  order  to  obtain  a  divergence 
free  result. 

In  expanded  form,  the  two  corrections  associated  with  the  sliding 
boundary  condition  are  found  to  be  as  follows? 

i  <k> 

■  if  (fi;)  {Aui:i(i,2,k)  +  Aun(i,2,k-1)  -  AuII(i-l,2,k)^uII(i-l,2,k-l) 

•Av11  (  "  '  -Avn(  "  )  -Avn(  "  )  -Av11  (  "  ) 

^w11  (  "  )-AwIX(  "  )-AwH(  "  )HAwn(  "  )}(B2-10) 

+  ■  i  ^2{  l^d^k)  +  — - 

-  ^  [^I(i-l,3,k)-»t!^p:i(l,3,k-l)-hCi<pI(i^l,3,ic)-»^pI(l,3,k+l)] 

■«AopI(i+l,3,k+l]  }  (B2-11) 

The  corresponding  results  for  the  upper  boundary  may  now  be  found 
by  a  systematic  permutation  of  superscripts  and  indices.  Hence  we 
obtain 
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(1) 

■  §L  1  [«^(i,N+2,k)+— ]  +  [«f£(i,N+3A)  +  — ]  } 

(2)  (3) 

-  (jjj)  {Aq>X(i,N+3,*)  +  — }  -  V2cpII(i,N+2>k) 

(4)  (5) 

•  -  (fiO  {  [f£(1,N+2,iO  +  —  ]  +  [F^(i,N+3>k)  +  —  ]  } 

(6)  (7) 

+  ~  Dn(i,N+2,k)  -  i  (^)  ^(1^+3^)  +  —  }  (B2-12) 


where  the  required  correction  terms  are 

(7) 

-  "  (fe)  {^(1^+3^)  +  —  } 

-  -  ^(Jl)  {AuI(i,N+3,k)-^uI(i,N+3,k+l)-AuI(i+l>N+3,k)-AuI(i+l,N+3,k+l) 

-  AvX(  "  )^vI(  "  )*vH  ”  )-AvI(  ”  ) 

+  Av!(  "  )*v1(  "  )-AvX(  "  )4AwI(  "  )}(B0-13) 

and  (2)  2 

“(fir)  {AflpJ(i,N+3,k)  +  — }  -  +  |  (£)  {  +  |AcpiI(l,N+2,k)  +  — 

+  [AcDII(i+l>N+2,k)+  A®II(i ,N+2,k+l)  +  AtpIX(i-l,N+2,k) 

+  A<pXX(i-l,N+2,k-l)]  -  (AcpIX(i+l,N+2,k+l)  +  A©XX(i+l,N+2,k-l) 

+  A(pII(i-l,N+2,k+l)  +  AtpXX(i-l,N+2,k-l)l  }  (B2-lh) 

This  completes  the  required  divergence  corrections  at  the  sliding 
boundaries. 
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APPENDIX  C 


EVALUATION  OF  KINETIC  ENERGY,  REYNOLDS  STRESSES,  VORTICITY  AND  SCALE 

OF  TURBULENCE 

C.l  REYNOIDS  STRESSES  AND  ENERGY 


In  differencing  method  C,  the  distributions  over  the  volume  element 
of  typical  velocity  components,  say  u  and  v,  are  assumed  to  be  of  the 
form  below,  for  control  volumes  of  either  family. 

u  =  AQ  +  A]_?  +  A2ti  +  A3C  +  A4T)C  +  A5C?  +  A6?t)  +  A7?tiC 

+  A8(i-?2)(i-n2)(i-f.2)  (ci-i) 

v  =  BQ  +  Bj_?  +  B2T1  +  B3C  +  B^ttC  +  B5C?  +  Bg^tl  +  B^ttC 

+  B8(1-?2)(1-t12)(1-C2) 

The  typical  mean  quadratic  product  uv  over  the  volume  element  is  then 
found  by  evaluation  of  the  integral 

+1  +1  +1 

uv  "  B  J  I  I  uv  d?dT^C  (Cl-2) 

-1  -1  -1 

Since  components  u  and  v  each  have  nine  terms  per  Eq.  (Cl-l),  the 
integral  in  (Cl-2)  has  8l  terms  in  all.  Fortunately,  most  of  these 
vanish  in  the  integration.  The  following  basic  types  of  integrals  are 
encountered  in  the  above  calculation,  namely 

1*5 


for  n  odd 


+1 


5  d?  =  0 


n+1 


for  n  even 


f  (i-t2)  a?  .  | 

-1 

+1 

f  (l-^2)  ?d?  =  0 

-1 

i  (i-?2)2  «  -  if 

-i 


(C1-3) 


and  similarly  for  the  variables  l"|  and  C» 

Consequently,  the  required  mean  product  reduces  finally  to  the 

form 

37  ’  {  AoB0  +  J  (fllBl  *  W  A3B3)  *  5  +  A5b5  +  A6b6> 

+  ^  (*7^7)  +  (a0b8  +  ;'0Bo)  +  3375  *8%  1  (cl-'*) 

The  first  term  on  the  right  represents  the  product  of  the  mean 
centroidal  values.  The  remaining  terms  reflect  the  fact  that  the  mean 
product  over  the  volume  differs  appreciably  from  the  simple  product  of 
the  two  centroidal  mean  values. 

The  overall  mean  product  is  then  found  by  averaging  of  the  above 
quantity  over  all  the  volume  elements  that  constitute  the  block.  More' 
over,  the  averaging  must  be  carried  out  over  the  volume  elements  of 
both  families.  Hence  the  total  volume  is  covered  twice.  This  is  a 
peculiarity  of  differencing  method  C. 


1U6 


By  means  of  the  usual  met  tod  of  permutation  ^  „ 

may  be  extended  to  include  all  six  of  th  '  6  reg°in«  re8ult 

ducts,  namely,  m  w  w - — 6  508611316  quadratic  mean  pro- 

stresses.  Furthermore,  the  sum  of  the  fir 
above  quantities  represents  twice  the  mean  k„et  -  ** 

this  summation  yields  the  result  ‘“"S'  ^  ln  f8ct 

*  ■  H  *  4  ♦  ]  ♦  i  [(**  ♦  ^  t  (B*  .  E=  +  b2} 


+  K  *  4  *  4,  ] 


’  5  f(fl2  +  4  *  +  (4  *  4  *  Bjj)  t  (Cjf  ♦  C|  ♦  c2,  ] 

+  57  l\  +  4  +  °7  I  +  H  I  Vs  +  B0BS  *  C0C8  1 
*  3#5 ' 4  *  DS  *  C§  | 


(ci-5) 
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C.2  VORTICITY 


\ 


The  three  vorticity  components  at  an  arbitrary  point  within  a 
given  control  volume  are  expressed  by  the  relations 


(C2-1) 


These  expressions  would  be  exact  if  the  derivatives  appearing  in 
them  were  known  exactly.  Unfortunately,  only  the  approximate  informa¬ 
tion  represented  by  equations  of  the  form  (Cl-l)  is  available.  Direct 
analytical  differentiation  of  (Cl-l)  degrades  the  accuracy  somewhat, 
although  the  subsequent  volume  integration  explained  below  tends  to 
offset  this  error.  The  results  obtained  are  of  the  form 

<is>  ®x  ■  <Vb3>  1  <C6-V«  •  V  -  V 

-  yn  +  yc  -  2c8(i-pi)ri(i -r2) 

+  2Bg  (l-?2)(l-n2)f,  (C2-2) 

The  components  and  uuz  follow  from  this  by  cyclic  permutation. 

The  mean  square  vorticity  finally  required  is  defined  by  the 
integral 

m  5  f  I  I  fax  +  ®y  +  *z  ]  drdnJ^  (C2_3) 

-1  -1  -1 
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Fortunately  most  of  the  integrals  involved  in  the  evaluation  of  (C2-3) 
vanish  identically.  The  integration,  although  quite  lengthy,  is  not 
difficult.  The  result  finally  obtained  is 

*  =  ^A3_CV  +  ^Bl"A2^  +  <W  1  +  3 

+  (B5-Au)c  +  (C6-35;2  +  A2  +  A26  +  Bb  +  e(  +  C2  +  C2  1 


+  I  ^  +  A  +  c?  1  +  ft!  ^8  +  4  +  C 


16 


+  27  ^a6B8  +  b6A8^  +  ^b1c8  ^  cUb8^  +  ^C5A8  +  A5C8^ 


This  result  must,  of  course,  be  averaged  over  all  the  volume 
elements  of  both  families,  in  the  same  manner  as  for  the  kinetic 
energy. 


C.3  SURFACE  CONSTANTS 


Consider  a  typical  control  volume  with  a  centroidal  point  of 
family  I.  The  eight  corners  and  the  centroid  may  be  designated  as 
follows : 


UPC 

un(i,J,k) 

up/  - 

uII(i-l,J-l,k-l) 

UA  ’ 

un(i-l,J,k) 

V  c 

uXI(i, J-l,k-l) 

“b- 

un(i,J-l,k) 

V  b 

uXI(i-l,J,k-l) 

uc" 

uH(i,  j,k-l) 

V  " 

uII(i-l,J-l,k) 

u0  = 

u1  (i,J,k) 

(C3-D 


Similarly,  consider  a  typical  control  volume  with  a  centroidal 
point  of  family  II.  The  nine  corresponding  points  may  be  designated 
as  follows 


up  -  uX(i+l,J+l,k+l) 

V  " 

uA  -  uX(i,,)+l,k+l) 

V  " 

Ug  -  uX(i+l, J,k+1) 

V  - 

uc  ‘  uX(i+l, j+l,k) 

V" 

U0  = 

«n  U,J>k) 

^(i^k) 

uZ(i+l,J,k) 

uI(i, J+l,k)  (C3-2) 

uJ(i> J,k+1) 


Now  in  either  of  the  above  cases  we  may  apply  Eq.  (Cl-l)  to  all  nine 
characteristic  points  of  the  control  volume,  thus  obtaining  nine  simul¬ 
taneous  equations.  These  are  easily  solved  for  the  nine  surface 
constants  Aq  through  Ag,  with  the  following  results. 


A0  =  5  {  +  UP  +  UA  +  u3  +  Uc  +  V  +  V  +  UA  ;  +  V  } 

Ai  =  5  {  +  UP  •  UA  +  “b  +  uc  ■  uc '  '  V  +  V  •  V  } 

A2  =  B  {  +  UP  +  ua  '  “b  +  uc  "  V  +  V  -  V  ■  V  } 

a3  e  H  {  +  UP  +  UA  +  “b  "  uc  +  V  "  V  ■  V  •  V  } 

A4  “  B  {  +  Up  +  UA  "  UB  '  Uc  '  V  •  V  +  UA  '  +  V  } 

A5  -  5  {  +  up  •  UA  +  “B  '  Uc  •  V  +  V  '  V  +  V  } 

A6  =  H  +  UP  -  UA  -  UB  +  Uc  +  V  ■  UB/  ■  V  +  V  } 

A7  =  5  {  +  UP  '  UA  ■  “B  ’  Uc  +  V  +  V  +  UA  '  -  V  } 

A8  “  u0  "  A0 

Similarly;  the  B's  are  computed  from  the  v's  and  the  C's  from  the  w's. 
Thereupon  the  mean  energy  is  calculated  from  (Cl-5)  and  the  mean  square 
vorticity  from  (C2-4). 
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C.4  SCAI £  OP  TURBUI£NCE 


\ 


Consider  a  hypothetical  velocity  perturbation  of  the  form 

u  «  P  sin  sin 

v  «  Q  sin  — sin  — j—  (c4-l) 

»  .  R  aln  ^  .in 

where  P,  Q,  R  and  l  are  arbitrary  constants. 

The  mean  kinetic  energy  E  is  then  given  by 

III 

2E  -  J  J  |  [u2+v2+v2]  dxdydz  -  J-(F^-KJ2+R2)  (cU-2) 

4  0  0  0 

The  vorticity  is  defined  by 

«.-(***) 

and  similarly  for  U)  and  (u  by  cyclic  permutation. 

y  z 

Upon  differentiating,  squaring  and  integrating,  we  obtain  for 
the  mean  square  vorticity 

—  Ill 

^  f  ‘  1  I  \  +  «L  +  \  "I  dxdydz  =  Q)2  |(P2-K}2+R2)  (CU-U) 

0  0“  6 

The  scale  of  turbulence  k  is  defined  by  the  equation 

^  -  oSK  «*-5> 

*  U) 
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where  C  is  a  constant  whose  value  we  wish  to  establish.  We  choose 
this  constant  in  such  a  way  that  the  hypothetical  perturbation  of 
(C^-l)  shall  be  deemed  to  have  a  scale  of  turbulence  X  which  is  equal 
to  its  wavelength  -t.  The  logic  of  this  choice  .  a  evident  from  inspec¬ 
tion  of  (C4-1).  Hence  we  substitute  from  (C^-2)  and  (C4-4)  into  (04-5) 
and  set  X  ■  l. 

The  result  obtained  thereby  is 


(C4-6) 


Note  that  the  arbitrary  amplitude  parameters  P,  Q,  R  cancel  from  this 
result,  so  that  the  constant  of  proportionality  given  above  is  indepen¬ 
dent  of  these  quantities. 

Equation  (c4-6)  is  now  taken  as  the  definition  of  the  scale  of 
turbulence  X,  and  is  regarded  as  applicable  to  all  types  of  perturba¬ 
tions.  It  may  be  interpreted  as  a  kind  of  root  mean  square  wavelength. 
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APPENDIX  D 


FOURIER  ANALYSIS  OF  TURBUIENCE  IN  A  UNIFORM  SHEAR  FLOW 


D.l  COORDINATE  SYSTEMS 


In  this  appendix,  all  quantities  are  non-dimensional  in  the  sense 
explained  elsewhere  in  this  report. 

Consider  an  arbitrary  point  which  moves  along  with  the  mean  flow. 
Let  x,  y,  2,  t  represent  the  ordinary  ti  ne-coordinates  of  the  point  in 
a  rigid  Cartesian  axis  system;  we  term  these  the  Eulerian  coordinates. 
Let  ?j  f|,  Ct  t  represent  the  corresponding  Lagrangian  coordinates  of 
the  point.  The  Lagrangian  coordinates  are  defined  in  terms  of  the 
Eulerian  coordinates  as  follows : 

5  "  xt  =  o 

^  *  yt  =  o  (D1'1) 

c  =  zt  =  o 

T  ■  t 

Thus  at  time  t  =  0,  both  sets  of  coordinates  coincide. 

Note  that  for  a  point  moving  with  the  mean  flow,  the  Lagrangian 
space  coordinates  ?,  t),  C  are  constants,  whereas  the  Eulerian  space 
coordinates  x,  y,  z  are,  in  general,  functions  of  time.  Conversely, 
for  a  point  fixed  in  space  and  not  moving  with  the  mean  flow,  coordin¬ 
ates  x,  y,  z  are  constants,  whereas  ?,  T],  C  are  functions  of  time. 


In  the  present  case  of  a  uniform  parallel  shear  flow  whose 
(dimensionless)  shear  rate  Q  equals  unity,  we  see  that  for  a  point 
moving  with  the  mean  flow,  the  relations  simplify  to 


X=?+Tlt  -  5+  yt 
y  =  t]  =  constant 
z  ■=  C  ■  constant 

t  =  T 


(?  *  constant) 


(Dl-2) 


Therefore  it  suffices  for  this  case  tc  take  x,  y,  z,  t,  as  the  Eulerian 
coordinates,  and  ? .  y,  z,  t  as  the  Lagrangian  coordinates. 

Consider  a  seu  of  simple  cubic  grid  points  in  the  fixed  or  Eulerian 
system  of  reference.  Let  space  and  time  points  be  defined  as  follows: 


L  t* 

\  '  5  (1 


ff> 


‘  S  (J 


L 

2 


■k  '  s  <k  - 


i) 

2 


.  T  , 

t  *=  W  \m 

m  M  ' 


I  =  1,2,3,“““ “N 


J  -  1,2,3,“—  N 
k  -  1,2,3, - N 


m  ■  1,2,3, - M 


(Dl-3) 


Let  us  establish  a  corresponding  set  of  Lagrangian  grid  points, 
which  move  with  the  mean  flow.  At  time  t  *=  0,  however,  let  the  moving 
lagrangian  grid  coincide  with  the  fixed  Eulerian  grid .  Let  indices  a, 

J,  k,  m  represent  space  time  points  in  the  moving  system.  In  particular, 
for  such  a  moving  grid  point,  we  have 

x  (o,J,k,m)  *  |  (a  -  |>  +  |  (J  -  |  -  §)  |  (">  -  §)  (Dl-4) 


where  the  indices  a,  J,  k,  m  vary  over  the  ranges  previously  indicated. 
The  y,  z  and  t  coordinates  remain  as  given  in  (Dl-3)* 
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D.2  FOURIER  SERIES  APPROXIMATION  FOR  VELOCITY 
PERTURBATIONS 


It  has  been  shown  elsewhere  in  this  report  that  it  is  permissible 
for  purposes  of  approximation  to  treat  turbulence  as  if  it  were 
periodic  in  space,  provided  that  the  wavelength  is  made  sufficiently 
large  and  the  cell  size  sufficiently  small.  Corresponding  considera¬ 
tions  apply  also  to  the  time  dimension.  In  principle,  therefore,  the 
space  time  structure  of  turbulence  is  representable  by  a  four- dimen¬ 
sional  Fourier  series  with  a  very  large  number  of  terms.  The  adequacy 
of  such  a  representation  in  practice  depends  on  the  number  of  terms 
which  it  is  feasible  to  retain,  and  also  on  the  accuracy  and  precision 
with  which  the  Fourier  coefficients  can  be  computed  and  expressed. 

We  introduce  the  auxiliary  variable 

6  -  ^  (p?  +  qy  +  rz)  +  (st)  (D2-1) 

where  p) 

q)  -  1,  2,  3,  — -  N 

r) 

s  -  1,  2,  3,  — -  M 


Now  the  u  velocity  component  is  expressible  in  the  form  of  the 
Fourier  series 


(D2-2) 
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In  this  equation  the  velocity  component  u  is  normalized  with 
respect  to  the  root  mean  square  turbulent  velocity  as  represented  by 

where  E  is  the  mean  turbulent  kinetic  energy.  The  variables  sin  9  and 
cos  0  are  also  normalized  by  their  respective  root  mean  square  values, 
both  of  which  equal  1//2.  Such  normalizing  is  optional,  but  does  confer 
upon  the  Fourier  coefficients,  the  A's  and  D's,  certain  advantageous 
properties,  as  will  be  seen.  The  matrices  A(p,q,r,s)  and  D(p,q,r,s) 
are  four- dimensional  arrays,  each  array  consisting  of  N^M  constants. 

Specifically  at  the  moving  Lagrangian  grid  points  themselves, 

Eqs.  (D2-1)  and  (D2-2)  reduce  to  the  forms 

6  =  f  fa  *  ib  +  +  r^k '  +  &r  Ls^m  ■  ] 


N  N  N  M 


sin  _0 

Il//2) 


(D2-3) 


Corresponding  expressions  for  the  other  two  velocity  components 
can  be  written  by  permuting  the  foregoing  equations  according  to  the 
pattern 


/\  /\  A 

w^ - v  C  <  ■■■  -  B  F< - ] 


(D2-4) 


Certain  features  of  the  above  results  are  noteworthy.  First  of 
all,  it  may  be  seen  that  once  the  Fourier  coefficients  (as  well  as  the 
normalizing  mean  energy  E)  are  known,  the  turbulent  velocities  are 
fully  defined.  The  Fourier  coefficients  are  the  six  arrays  of  constants 
namely,  the  A's,  B's,  C's,  D's,  E's  and  F's.  The  total  number  of  Fourier 
coefficients  is  therefore  6n^M,  or  just  twice  the  total  number  of 
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velocity  components  in  the  complete  space  time  grid  block.  Theoretically, 
these  coefficients  define  the  turbulent  motion  to  Just  that  degree  of 
resolution  permitted  by  the  constants  N  and  M  which  fix  the  fineness  of 
the  coordinate  mesh.  The  adequacy  and  resolution  of  this  mathematical 
model  theoretically  improve  as  the  constants  N  and  M  are  increased. 
However,  as  the  number  of  Fourier  coefficients  increases,  the  number  of 
significant  figures  required  for  each  constant  also  increases,  in  order 
to  prevent  difficulties  arising  from  round-off  errors. 

Of  course,  the  velocity  components  defined  by  Eqs.  (D2-2)  or  (D2-3) 
are  referred  to  the  moving  Lograngian  coordinate  system.  However,  if  a 
transformation  to  fixed  Eulerian  coordinates  is  required,  this  is  easily 
found  by  use  rfEqs.  (Dl-2).  The  reason  for  introducing  Lagrangian 
coordinates  in  the  first  place  is  that  the  sliding  boundary  conditions 
along  the  top  and  bottom  of  the  rigid  Eulerian  block  are  thereby  satis¬ 
fied. 

It  should  be  stressed  that  the  above  Fourier  representation  is 
complete  in  the  sense  that  it  duplicates  the  actual  turbulence  to 
within  the  allotted  degree  of  resolution.  Hence  any  or  all  statistical 
or  phenomenological  quantities  that  can  be  calculated  from  a  finite- 
difference  approximation  of  the  velocity  history,  can  also  be  calcu¬ 
lated  directly  from  the  Fourier  coefficients  themselves. 

Inasmuch  as  the  Fourier  coefficients  depend  only  on  the  wave 
numbers  j?,  £  and  |  but  not  on  the  space  variables,  ?,  y,  z,  the 

turbulence  is  homogeneous.  Also,  since  the  Fourier  coefficients  are 
independent  of  the  time  variable  t,  the  turbulence  is  stationary. 

An  important  advantage  of  the  Fourier  representation,  as  compared 
with  the  direct  velocity  components  themselves,  is  that  it  more  clearly 
reveals  the  structure  of  the  turbulence.  The  direct  velocity  components 
in  physical  space  seem  so  chaotic  as  almost  to  defy  any  efforts  to 
analyze  or  order  them.  On  the  other  hand,  the  Fourier  spectral  coef- 
ficie  its  can  be  expected  to  vary  in  some  smooth  and  orderly  fashion 
within  the  wave  number  domain.  Hence  there  can  be  some  hope  of  charting 


158 


this  spectral  variation,  and  perhaps  of  expressing  it  graphically  or 
analytically  in  some  greatly  condensed  and  simplified  form.  If  this  can 
Ibe  accomplished  it  will  mean  that  the  problem  of  turbulence  In  a  uniform 
ehear  flow  is  then  essentially  solved.  However,  since  the  wave  number 
domain  is  four-dimensional,  to  summarize  conditions  in  this  realm  in  a 
way  which  is  both  simple  and  adequate  may  or  may  not  be  possible.  De¬ 
tailed  numerical  data  is  required  before  this  question  can  be  settled. 

Attention  is  invited  to  the  fact  that  the  Fourier  solution  contains 
equal  numbers  of  sine  and  cosine  terms.  From  this  fact,  it  follows  that 
the  overall  turbulence  i6  resolvable  into  components  at  each  point  of 
wave  number  space  which  are  in  certain  sense  out  of  phase  or  orthogonal 
to  each  other.  This  is  important  for  the  adequate  description  of  turbu¬ 
lence.  Conventional  experimental  methods  of  obtaining  spectral  informa¬ 
tion  about  turbulence  usually  fail  to  disclose  these  phase  relationships, 
and  are  therefore  incomplete  in  a  very  significant  respect,  (in  fact  for 
a  complete  description,  information  on  phase  relationships  is  Just  as 
essential  as  information  about  amplitudes;  in  quantitative  terms  we  may 
assert  that  exactly  half  the  needed  information  concerns  amplitude,  the 
other  half  concerns  phase.) 

In  principle,  the  equation  of  continuity  and  the  three  equations  of 
motion  could  be  written  directly  in  terms  of  the  Fourier  coefficients 
themselves  rather  than  in  the  more  usual  way,  in  terms  of  velocity  compon¬ 
ents.  These  equations  could  then  be  solved  for  the  'unknown  Fourier  coef¬ 
ficients.  This  method,  however,  seems  to  involve  extreme  analytical  com¬ 
plexities  and  numerical  difficulties,  and  has  been  discarded  as  being  too 
unpromising.  On  the  other  hand,  a  feasible  method  of  solution  has  been 
developed  which  establishes  the  velocities  directly,  without  any  reference 
to  Fourier  coefficients.  Thus  far  in  this  discussion  we  have  considered 
how  unknown  velocities  may  be  expressed  in  terras  of  known  Fourier  coeffici¬ 
ents.  In  the  present  circumstances,  however,  the  velocities  are  the  knowns 
and  the  Fourier  coefficients  the  unknowns.  Hence  we  face  the  problem  of 
inverting  the  previous  solution.  Fortunately,  it  is  characteristic  of 
Fourier  methods  that  the  fundamental  equations  can  be  readily  inverted. 
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D.3  SOLUTION  FOR  THE  FOURIER  COEFFICIENTS 


\ 


For  the  purpose  of  inverting  the  basic  equation  (D2-2),  we 
introduce  the  following  notation 


0  -  (p^  +  qy  +  rz)  +  y  (st) 


8'er(p/?M/y  +  r'z)  +  =r  (s't) 


2it 


T  L  L  L 


_  f  f  f  f  sin  9  sin  0;  d?  dy  dz  dt 
“  \io  J  yio  4  ^  ^  L3T 


T  L  L  L 


(D3-1) 


f  !*  f  I  cos  6  cos  0;  d?  dy  dz  dt 

XC=  .0,  .1  .i  tI77f)TWT  ~3t' 


t=0  z=0  y=0  ?=0 

T  L  L  L 


T  f  f  f  f  sin  9  cos  9 '  d%  dy  dz  dt 

zbo  “  J  J  J  -J  TU7T)  (i /JZ)  l3t 

t=0  z«=0  y=0  5=0  L  i 


It  can  be  verified  by  direct  integration  that  the  following  ortho¬ 
gonality  relations  are  valid,  namely, 


(a) 


ss 


D  -  p 

+1  provided  )  q  *  «=  q 
•'  -  r 


s  «=  s 


0  othervise 
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(b) 
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p 
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I  =  +  1  provided 


/ 

P  ■ 


r  =  r 


s  =  s 


(D3-2) 


=  0  otherwise 

(c)  I  =0  for  all  values  of  the  indices. 


sc 


We  now  multiply  Eq.  (D2-2)  through  by  sln  ?  dL dy_d 2  gt 

(1//2)  L3T 


and  integrate  over  the  limits  indicated  in  (D3-1). 
orthogonality  relations  (D3-2),  we  obtain 


Upon  invoking  the 


I  J 


t=0  z=0  y=  0  5=0 


■  z,t)  sin  6  d5  dy  dz  dt 

f2E  (l//2)  L'T 


A(p,q,r,s)  (D3-3) 


For  purposes  of  approximate  numerical  evaluation,  the  integrals  on  the 
left  may  be  replaced  by  corresponding  finite  summations.  These  are 
fonned  by  multiplication  of  each  principal  grid  value  of  u  by  the  cell 
volume  (jO3  and  time  interval  (£) ,  and  then  summation  over  the  full  extent 
of  the  space  time  block.  Upon  rearranging,  we  obtain  the  result 


N  N  N  M 


A(p,q, 


■  htS  5  S  5"“$ 


,m)  sin  8 

(1//2) 


0>3-M 


Proceeding  in  like  fashion  with  the  cosine  term  we  obtain  the  corres¬ 
ponding  result 


-s  i  it 


(03-5) 


l6i 


The  remaining  Fourier  coefficients  may  be  found  in  like  manner  by  permu¬ 
tation  in  accordance  with  (D2-4). 

The  velocity  components  on  the  right  side  of  Eqs.  (D3-*0  and 
(D3-5)  are  evaluated  at  the  Lagrangian  grid  points.  These  nu/  be 
obtained  by  interpolation  from  the  previously  computed  and  known  values 
at  the  Eulerian  grid  points. 

A  further  important  result  is  obtained  from  energy  considerations. 
First,  both  sides  of  Eq.  (D2-2)  are  squared  and  integrated  over  the 
space  time  block.  In  this  process,  the  right  side  is  greatly  simpli¬ 
fied  through  the  application  of  the  orthogonality  relations.  The  same 
procedure  is  also  applied  to  the  other  two  components,  and  the  three 
equations  are  then  added  together.  The  following  result  is  thereby 
obtained,  namely, 

Y.  Y  Y  S  {  |"A2(p»q>r'8)  +  D2(P><3,r,s)] 

p»l  q*l  r»l  6=1 

+  [B2(p,q,r,s)  +  E2(p,q,r,s)  ]  +  [c2(p,q,r,s)  +  F^(p,q,r,s)j  }  =  1  (D3-6) 

The  physical  interpretation  is  simple.  Each  quadratic  term  represents 
the  fraction  of  the  total  kinetic  energy  of  turbulence  stored  at  the 
corresponding  location  in  wave  number  space.  The  sum  of  all  such 
eneigy  fractions  must  equal  unity.  This  simple  form  of  result  is  the 
consequence  of  the  normalizing  procedures  employed  in  the  analysis. 

It  is  perhaps  worth  mentioning  that  a  shift  in  the  origin  of  the 
Lagrangian  space  time  coordinates  will  affect  all  of  the  Fourier 
coefficients,  but  in  a  simple  way.  The  total  energy  stored  at  each 
point  in  wave  number  space  remains  unchanged,  but  the  separate  sine 
and  cosine  components  each  undergo  a  simple  phase  shift. 

It  might  also  be  mentioned  that  all  results  obtained  in  this 
appendix  by  use  of  Fourier  series  have  an  equivalent  representation  in 
terms  of  truncated  Fourier  integrals.  The  difference  is  mainly  one  of 
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form,  the  series  form  being  somewhat  more  convenient  for  numerical 
methods.  In  the  integral  method,  the  Fourier  coefficients  are  treated 
not  as  discrete  constants  but  as  continuous  density  functions  in  wave 
number  space.  Conversely,  the  discrete  Fourier  constants  may  be 
interpreted  as  pointwise  approximations  to  functions  which  are,  in 
fact,  distributed  continuously  in  wave  number  space.  The  distinction 
is  largely  academic.  There  is  no  difference  in  the  final  computational 
formulas  which  are  obtained ! 

This  completes  the  essential  features  of  the  Fourier  analysis.  It 
is  now  possible  to  derive  furtner  interesting  and  useful  results  con¬ 
cerning  quantities  such  as  the  Reynolds  stresses,  various  correlation 
coefficients,  and  so  on,  but  such  additional  topics  lie  beyond  the 
scope  of  the  present  discussion. 

As  a  final  cautionary  note,  we  again  observe  that  calculation  of 
the  Fourier  coefficients  from  Eqs.  (D3-^)  and  D3-5)  is  subject  to 
inaccuracy  from  round-off  errors,  especially  for  high  values  of  N^M. 

It  will  be  necessary  to  retain  a  sufficiently  large  number  of  signifi¬ 
cant  figures  in  the  calculations.  It  is  a  curious  fact  that  if  N^M 
is  small,  round-off  error  in  the  individual  coefficients  is  small, 
but  detailed  resolution  of  the  turbulence  is  poor.  On  the  ether  hand 
if  N^M  is  high,  resolution  is  good,  but  round-off  error  is  high.  This 
circumstance  somewhat  resembles  the  well  known  indeterminacy  principle 
of  Heisenberg! 
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never  before  possible.  In  order  to  accomplish  this  aim,  however,  it  is  firtls,  nec¬ 
essary  to  formulate  the  essential  theoretical  concepts  in  a  suitable  manner.  This 
report  summarizes  the  progress  achieved  to  date  in  this  connection.  Various  es¬ 
sential  basic  equations  are  derived,  but  the  emphasis  is  as  much  on  fundamental 
concepts  as  on  mathematical  details.  More  specifically,  a  method  is  established 
for  the  computer  simulation  of  the  detailed  stationary  turbulence  in  a  uniform 
shear  flow.  The  results  obtainable  in  this  way  are  far  more  comprehensive  than 
any  which  could  reasonable  be  obtained  by  physical  experiment.  The  data  generated 
represent*  fundamental  information  which  may  be  subsequently  analyzed  to  establish 
overall  phenomenological  characteristics  of  the  turbulence.  The  concepts  in  this 
report  should  provide  a  sound  basis  for  a  systematic,  sustained  and  productive 
research  plan.  They  have  already  been  successfully  applied  to  a  computer  program 
which  is  now  going  into  operation.  Results  of  a  typical  computer  run  are  included 
and  illustrate  qualitative  agreement  with  theoretical  predictions;,  it  is  hoped 
to  present  far  more  comprehensive  and  definitive  numerical  results  in  future  re-  | 
port 8. 
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